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and all influence coefficients are derived in closed form. These and other features combine to 
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realities of a user-oriented environment. A wide variety of numerical results demonstrating 
the method are presented. 
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1.0 SUMMARY 


A general panel method for solving the boundary value problems of linear subsonic potential 
flow in three dimensions is presented. The method is characterized by a building block 
approach wherein an influence coefficient equation set is developed by assembling panel 
networks appropriate to the boundary value problem. Each network is composed of a 
paneled portion of the boundary surface on which a source or doublet singularity distribution 
is defined, accompanied by a properly posed set of boundary conditions. Curved panels 
possessing linearly varying source or quadratically varying doublet singularities are employed, 
and all influence coefficients are calculated in closed form. Both analysis (Neumann) and 
design (Dirichlet) boundary conditions are treated. In this treatment standard aerodynamic 
auxiliary conditions (e.g., Kutta, closure, and continuity conditions) arise as natural boundary 
conditions associated with various network types. 

A pilot computer program was developed to assess the feasibility of the method and a wide 
variety of check cases were run. These included thin and thick wing analysis and design 
problems, wing-body analysis and design problems, convergence studies, trade studies, and 
timing checks. These check cases indicate that the method has the necessary flexibility, sta- 
bility and accuracy, and efficiency. 

1. Flexibility: The method offers the user a variety of modeling options including actual 
and mean surface paneling as well as velocity and potential boundary conditions. 

2. Stability and Accuracy: The method displays a marked insensitivity to the size, shape 
and arrangement of panels and achieves accurate results with relatively sparse panel 
densities. 

3. Efficiency: The method appears to be almost as efficient as existing first order 
techniques on a panel by panel basis. On acase-by-case basis the present method has 
significant advantages because it requires fewer panels than first order methods. These 
advantages become extremely important for complex configurations. 



2.0 INTRODUCTION 


Steady, inviscid, irrotational, and incompressible fluid flow in a domain D is characterized by 
a perturbation velocity potential <j>(P) satisfying Laplace’s equation 


•^xx + <^yy + ^zz ■ 0 

(compressibility effects over a wide range of subsonic Mach numbers can be approximated 
by the Gothert rule in which case a coordinate transformation again results in Laplace’s 
equation ). 

Fluid flow boundary conditions associated with Laplace’s equation are generally of analysis 
or design type. Analysis conditions are employed on portions of the boundary where 
geometry is considered fixed, and resultant pressure distributions are desired. The permea- 
bUity of the fixed geometry is known; hence, analysis conditions are of the Neumann type 
(specification of normal velocity). Design boundary conditions are used wherever a geometry 
perturbation is allowed for the purpose of achieving a specific pressure distribution. Here a 
perturbation to an existing tangential velocity vector field is made; hence, design conditions 
are fundamentally of the Dirichlet type (specification of potential). The design problem in 
addition involves such aspects as stream surface lofting (i.e., integration of streamlines 
passing through a given curve), and the relationship between a velocity field and its potential. 

Green’s theorem (ref. 1 ) shows that any solution of Laplace’s equation can be expressed as 
the potential induced by source and/or doublet singularities distributed on the boundary B 
of D with the singularity strengths being determined by the boundary conditions. This fact 
has been exploited for some time in the design of numerical solution techniques for use on 
large-scale digital computers. 

The Douglas Neumann program (ref. 2) was spectacularly successful for its time in solving 
complex potential flow problems with Neumann boundary conditions. The numerical method 
represented the boundary by constant strength source panels with the strengths determined 
by an influence coefficient equation set relating the velocities induced by the source panels 
to the boundary conditions. The lack of provision for doublet panels limited the class of 
solutions to those without potential jumps and hence without lift. 

On the otherhand, lifting solutions had long been generated by vortex lattice techniques using, 
in effect, constant strength doublet panels (ref. 3). These schemes were developed primarily 
for the analysis and design of thin wings. The lack of provision for source panels made the 
treatment of thick configurations difficult. 

One of the first computer program systems for attacking arbitrary potential flow problems 
with Neumann boundary conditions (refs. 4, 5, and 6) combined the source panel scheme 
of the Douglas Neumann program with variations of the vortex lattice technique to form a 
general boundary value problem solver (known as the Boeing TEA 230 program). In addition 
to this method, many other general schemes have been developed, typical of which are the 
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methods of references 7,8,9, 10, 11, and 1 2. These methods are all similar in technique, 
differing primarily in the specific integral equation used and/or the construction of the singularity 
distributions employed in its solution. Each scheme has achieved considerable success in a 
variety of applications. 

A very useful feature of the TEA 230 program is the ability to handle, in a logical fashion, 
any well-posed Neumann boundary value problem. From its inception, the method has 
employed a building block approach wherein the influence coefficient equation set for a 
complex problem is constructed by simply assembling networks appropriate to the boundary 
value problem. A network is viewed as a paneled surface segment on which a source or 
doublet distribution is defined, accompanied by a properly posed set of Neumann boundary 
conditions. The surface segment can be oriented arbitrarily in space and the boundary 
conditions can be exact or linearized. Several doublet network types with differing singularity 
degrees of freedom are available to simulate a variety of physical phenomena producing 
discontinuities in potential. These features combine to allow the analysis of configurations 
having thin or thick wings, bodies, nacelles, empennage, flaps, wakes, efflux tubes, barriers, 
free surfaces, interior ducts, fans, etc. 

While the extreme versatility of TEA 230 has been well appreciated over the past decade, 
the need for basic improvements has become clearly evident. Some of these stem from the 
fact that it sometimes takes weeks requiring the expertise of an engineer having years of 
experience with the method to set up and run a complex configuration. To some extent 
this is unavoidable; in order to correctly model a complex flow for which no prior user 
experience is available, the engineer must understand the properties and limitations of 
potential flow. However, once the boundary value problem has been formulated, the user 
must still contend with certain numerical idiosyncrasies and inefficiencies which require 
adherence to stringent paneling rules-rules which are frequently incompatible with the 
complex geometrical contours and rapidly changing aerodynamic length scales of the vehicle 
under analysis. These difficulties are directly related to the use of flat panels with constant 
source and doublet strengths. Methods employing these features seem to be quite sensitive 
to panel layout. Numerical problems arise when panel shapes and sizes vary, and fine 
paneling in regions of rapid flow variations often forces fine paneling elsewhere. In addition, 
large numbers of panels are required since numerical accuracy is often strongly affected by 
local curvature and singularity strength gradient. Such problems place severe limitations on 
the development of automatic panelers or other complementary aids aimed at relieving the 
user of the large amount of handwork and judgements associated with producing accurate 
numerical solutions involving complex geometrical shapes. 

The versatility of the Boeing TEA 230 program in a user oriented environment has 
motivated the adoption of a similar, but less sensitive building block approach for the present 
method. In fact, the same network concept has been adopted and generalized to include 
Dirichlet boundary conditions. The treatment of Dirichlet boundary conditions not only 
provides the capability for designing surface segments to achieve desired pressure distributions, 
but also clarifies the nature of the boundary value problem associated with modeling viscous 
wakes and other effects through the introduction of discontinuities in the potential. 
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In addition, the present method has sought to enhance practical usability by improving upon 
the flat, constant singularity strength panels employed in the construction of the networks. 
The details of the present method arose from a prior analysis identifying the numerical 
features required to eliminate the practical difficulties encountered by the TEA 230 program. 
These features include the use of curved panels, linearly varying source strengths and 
quadratically varying doublet strengths. The present method implements these features in 
such a manner that all influence coefficients are evaluated in closed form thereby avoiding 
numerical integration— a weak link where singular integrals are involved. Other higher order 
panel methods have recently been studied and in some cases implemented (refs. 9, 13, and 14) 
with excellent results, but the present method is the first to cover in comprehensive fashion 
the complete spectrum of possible analysis and design boundary value problems. 

The author wishes to acknowledge that section G. 1 and portions of appendix D were pre- 
pared by Larry L. Erickson, NASA-Ames Research Center. 
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3.0 ABBREVIATIONS AND SYMBOLS 


a 

a 

^j 

JR 

A 

b 

B 

B 

B' 

B" 

c 

'c 


c(t?) 


Cp 

CpY 

^FZ 


defined by equation (A. 20) 

coordinate in direction of v defined in figure D.3 (meters) 
elements of A defined by equation (D. 149) 
wing aspect ratio 

matrix of influence coefficients; also, orthogonal coordinate transformation 
defined by equation (A.5) 

vector of specified boundary conditions; also, defined by equation (A. 20); 
also, wing span (meters) 

element of b 

boundary of D 

matrix relating panel source or doublet distribution coefficients to singu- 
larity parameters in a neighborhood of the panel 

portion of B 

portion of B 

defined by equation (D.24); also, wing chord (meters) 

defined by equation (C7) 

defined by equation (C.7) 

defined by equation (C.3) 

defined by equation (C.3) 

reference chord (meters); also, column matrix of source or doublet 
coefficients 

section chord (meters) 

drag coefficient 

force coefficient 

force coefficient in Y direction 

force coefficient in Z direction 

lift coefficient 
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<L 

C(M, N) 

CMx 

^My 

CP 

Cp 

‘'Po 

CR 

C], C2 

d 

dp 

‘^H 

dp (^2 

D 

D(M, N) 

D' 

E(M, N, K) 
E1,E2, E3,E4 
F(M, N, K) 
g 

A 

g 

G(M, N) 
h 

H(M,N, K) 


centerline 

panel far field moments defined by equation (D.83) 

moment coefficient about X axis 

moment coefficient about Y axis 

central processor 

pressure coefficient 

existing or initial pressure coefficient 

reference chord (meters) 

defined by equation (G.18) 

tail span (meters) 

defined preceding equation (D.60) 

defined preceding equation (D.41) 

defined by equation (G.22) 

fluid domain 

defined by equation (D.98) 
domain adjacent to D 

panel edge endpoint functions defined by equation (D.58) 

recombinations of far field moments 

panel edge line integrals defined by equation (D.40) 

/iW 

see equation (D.141) 

defined by equation (D.95) 

defined by equation (D.15) 

panel surface integrals defined by equation (D.25) 
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H*(M,N,K) 


KM, N) 
J(M, N) 

C 

A 

C 

Ci,82 

L 

M 

MXK 

MXQ 

MXFK 

My 


N 

NHK 


P 

P 



complementary H integrals defined by equation (D.51 ) 
defined by equations (D.23), (D.32), (D.77) or (D.87) 
defined by equation (D.29), (D.34), (D.85) or (D.92) 
defined by equation (C. 1 1 ) 
defined by equation (C.l 1 ) 

panel edge coordinate defined in figure D.3 (meters) 
unit tangent vector along edge of 2 
values of £ at edge endpoints 
typical side of S 

number of comer point rows (appendix A) 
defined by equation (D.36) 
defined by equation (D.36) 
defined by equation (D.54) 

moment coefficient about axis parallel to Y axis through reference point 
freestream Mach number 
number of panels 
unit surface normal 

number of comer point columns (appendix A) 
defined by equation (D.49) 

X coordinate of center of pressure (meters) 
field point (meters); also, magnitude of ^ (meters) 
field point position vector 
centroid of 2 

panel center defined by equation (A. 2) 


Pj , P2, P3, P4 panel comer points 





Q 

-► 

Q 

r 

R 

R 

Rr 

S 

Sr 

si,S2 

♦ 

t 

-»■ 

-► 

tL 

tu 

A 

t 

T 

Tr 

V 
7 

V 

Vt 


p/p 

point on boundary B (meters); also, magnitude of also, point on panel S 
panel point position vector 
distance from P to Q (meters) 

weighted residual function; also, defined by equation (D.2) 
position vector (x, y, z) (meters) 
origin of local panel coordinate system 
origin of moment 

simply connected portion of B; also panel surface 

'y 

reference area (meters^) 
defined by equation (G.18) 
tangential vector field on S 
defined by equation (C8) 
defined by equation (C8) 
defined by equation (C.4) 
defined by equation (C.4) 
unit vector in direction T 
superscript denoting transpose 

reference length for computing moments about axis (meters) 

defined by equation (D. 1 1 3 ) or (D. 1 3 1 ) 

perturbation velocity (meters/sec.) 

— ► 

magnitude of V (meters/sec) 
magnitude of (meters/sec) 
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V 


V 


A 





X 


(x, y, z) 


(xq, yo> zq^ 

y 

Z 


magnitude of (meters/sec) 
total fluid velocity vector (meters/sec) 

average of upper and lower surface velocities (meters/sec) 

defined by equation (D.l 14) or (D.132) 
defined by equation (D. 1 1 5) or (D. 1 33) 

defined by equation (D.134); also, difference at upper and lower surface 
velocities (meters/sec) 

lower surface velocity (meters/sec) 

tangential component of i.e., ^ -(It • (meters/sec) 

upper surface velocity (meters/sec) 

freestream velocity vector (meters/sec) 

chordwise coordinate (meters) 

field point coordinates in local panel coordinate system 
defined by equation (D. 1 5) 
spanwise coordinate (meters) 
vertical coordinate (meters) 


P 

Pt 

Pi’ Pi 

P<t> 

r 


angle of attack (degrees or radians) 
defined by equation (G.12) 
defined by equation (C.3) or (C.7) 
defined by equation (C.4) or (C.8) 
defined by equation (G.17) 
defined by equation (C. 1 1) 

set of control points; also, circulation (meters^/sec); also, gamma function 
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A 

ACp 

S 

5g 

6h 

e 


e(M, N, K) 

f 

V 

T?V 

6 

X 

X({, T7, f ) 

/ 

M 


jump in value across surface 
jump in Cp across surface 
defined by equation (A.21) 
defined preceding equation (D.60) 
defined preceding equation (D.41) 

defined by equation (D.17); also, defined by equation (D.72); also, 
spanwise gap (meters) 

does not belong to set 

defined by equation (D.52) 

defined by equation (A. 19) 

defined by equation (A. 19); also, tangential coordinate normal to panel 
edge (meters); also, fractional semispan 

fractional semispan (horizontal) 

fractional semispan (vertical) 

polar angle (degrees) 

vector of unknown singularity parameters 
element of X 

singularity parameter 
singularity distribution on a panel 
doublet strength (meters /sec) 

tangential derivation of doublet strength normal to a panel edge (meters/sec) 


(Mq» panel doublet distribution coefficients 

(Mq> My Mxx» Mxv> Myy) panel doublet distribution coefficients expanded about 

the point (x, y, 0) 

u (M, N, K) defined by equation (D.52) 

0 = (i^j, unit edge normal defined in figure D.3 

^ defined by equation (A. 1 9) 

(^, T 7 , f) panel point coordinates in local panel coordinate system 
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AAA 

(hv. n 






3/3np 


orthogonal unit vector defining local panel coordinate system 
defined by equation (D.19) 
source strength (meters/sec) 

panel source distribution coefficients expanded about the point (x, y, o) 
panel source distribution coefficients 

plane quadrilateral formed by projecting panel comer points onto the 
local (?, 7 ]) plane, which for the flat approximation is the same as the 
near plane 

summation over four panel edges 

polar coordinate defined by figure G.l ; also, perturbation velocity poten- 
tial (meters^/sec) 

average perturbation potential (meters^/sec) 
difference perturbation potential (meters^/sec) 
upper surface perturbation potential (meters /sec) 
lower surface perturbation potential (meters^/sec) 
total velocity potential (meters^'/sec) 
existing or initial velocity potential (meters /sec) 
angle of yaw (radians) 
set of singularity parameter locations 


leading edge sweep angle (degrees); also, set of all singularity parameters 


A ^ 
n • Vp 

A 

n . Vq 

desired change in Cp 
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V 

Vp 

Vq 

(^) 

D 


gradient operator 

gradient operator with respect to coordinates of P 
gradient operator with respect to coordinates of Q 
denotes a vector, e.g., the position vector P 


generally denotes a matrix, e.g., P = 
elements are the coordinates of P 
denotes vector cross product 



is the column matrix whose 
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4.0 THEORY 


In this section we develop the integral equations which are the basis of the solution technique 
of the present method. In section 4.1 we consider the fluid flow boundary value problems 
associated with equation (1) and then in section 4.2 we use Green’s theorem to derive 
equivalent integral equations. 

4.1 FLUID FLOW BOUNDARY CONDITIONS 

In this section we consider the specification of boundary values for the perturbation poten- 
tial (f> of equation (1 ). The boundary B of the fluid domain D contains of course the sur- 
face of the configuration being analyzed or designed. To account for viscous wakes, B is 
often augmented by artificial wake surfaces with boundary conditions which allow for po- 
tential jumps. In particular, a domain D which is originally multiply connected is usually 
made simply connected in this manner, thereby allowing the imposition of Kutta or other 
conditions (see fig. 1 ). The placement of these artificial surfaces and the selection of appro- 
priate boundary conditions is a physical modeling problem of some complexity. We do not 
address such a problem specifically as it is beyond the scope of this report. However, we do 
consider boundary conditions which are sufficiently general to allow for the treatment of a 
wide variety of physical models. 



{Note: For simplicity, we display a two rather than three-dimensional domain) 


Figure 1. — Fluid Domain and Boundary 
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As a general rule, a boundary value problem associated with equation (1) is well posed if 

00 

either <p or its normal derivative — is specified at every point of B (ref. 1 , chapter 12). 

0n 

(For simplicity, we do not address the regularity assumptions required for a rigorous treat- 
ment of the subject.) One major exception is the interior Neumann problem where D is a 

00 

finite domain and only — is specified on B (fig. 2). Gauss’ theorem shows that a solution 
can exist only if 


-ds = 0 (2) 

0n 

B 

00 

In case the specified values of — satisfy equation (2), an infinite number of solutions for 

0n 

0 exist, although they differ by an additive constant only. If D is an infinite domain, 

(fig. 1 ) the behavior of 0 near infinity must be restricted for a unique solution; however, 
sufficient restriction is an automatic consequence of the integral equation formulation of 
the next section, (see ref. 15, chapter 3). 




Figure 2, — Neumann Problem for Finite Domain 
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Neumann boundary conditions (specification of — ) arise naturally in the analysis of fixed 

0n 

configurations bounded by surfaces of known permeability. For example, let S be a por- 
tion of B bounding D on the side witl^surface normal n (fig. 3). If S is impermeable, the nor- 
mal component of the total velocity V must vanish on S. By definition, V is the gradient 
of the total potential ^ defined by 


<!> = 0 + R . 


(3) 


where R is the position vector (x,y,z) and is the freestream velocity vector. Hence, 

V=Voo+v0 (4) 

and the Neumann boundary condition expressing impermeability is 


where 


30 

3n 



30 , A 

— = v0 • n 
3n 


Here, V0 is the perturbation velocity vector. 


(5) 


( 6 ) 


n 
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Once a solution <l> to any boundary value problem has been found, a pressure coefficient 
Cp at each point on S may be computed by the formula 

Cp=l-(V/V«)2 (7 


where is the magnitude of and V is the magnitude of V . For the particular case 
of an impermeable boundary, equation (7) becomes 

Cp= 1 -(Vj/V„)2 (8) 


where is the magnitude of the tangential projection of V. 

Dirichlet boundary conditions (specification of <j>) arise in connection with the inverse prob- 
lem, that of solving for a specified pressure distribution on S. However, the achievement 
of a desired pressure distribution on S is not physically significant without restrictions on 
the flux through S. To achieve both specified pressure and normal flow distributions on S, 
the position of S must in general be perturbed. To elaborate further, we must assume (for 
simplicity) that S is required to be a stream surface. The specification of 4> guarantees a 
predetermined tangential velocity vector field, a-fortiori a predetermined pressure coefficient 
distribution on S via equation (8). However, the resultant normal velocity on S will not 
in general, satisfy equation (5) and for this, a modification of the normal vector distribution 
(hence position) of S is required. The total design problem is thus composed of two prob- 
lems. The first is to find a perturbation potential on S yielding a specified distribution of 
pressure coefficient as defined by equation (8), and the second is to update S to be a stream 
surface in the resultant flow. The two problems are coupled and in general, an iterative pro- 
cedure is required for solution. Sophisticated iteration techniques are available for this pur- 
pose (ref. 16). However, in most cases, the coupling is surprisingly weak and the two prob- 
lems may be solved sequentially in such a manner that (5) and (8) hold simultaneously after 
relatively few iterations. 

To solve the first problem, we assume that an approximation to the solution already exists. 
Such an approximation can be the result of a previous iteration or be defined initially by 
guessing a location of S and solving the pure analysis problem. The problem then, is to 
modify the existing potential so that its tangential velocity field yields a more desirable pres- 
sure coefficient distribution. For this purpose, we prefer to deal directly with the tangential 
velocity field itself. Our approach is to assume that 0 is continuous on S and to specify a 
vector field t on and tangent to S such that 

V • t = |t| where t = t / |t| (9) 

For simplicity, we assume that ^ is simply connected ^d that the field t is continuous and 
iTl ^ 0. (T is otherwise arbitrary). A good choice for t can be obtained by linearizing equa- 
tion (8) about existing values, i.e.. 
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where Cpo and are the existing values of pressure coefficient and tangential velocity 
respectively and 3Cp is the desired change in pressure coefficient. (We assume 
dCp/2( 1 - Cp^) is small compared to 1 .) 

Once a field t has been selected, certain constants of integration remain to be specified be- 
fore equation (9) defines actual Dirichlet boundary conditions. One possible choice is the 
specification of 4> (and hence 0) at all points of inflow of the fieldTinto S, i.e., points on the 
boundary of S where t is directed into S (see fig. 4). This choice allows the integration of 
on S and is thus directly equivalent to Dirichlet boundary conditions. The specified values 
of 4> should be close to the existing values 4 >q, otherwise the resultant tangential velocity 
will not be close to t and consequently Cp will not be close to Cpo + 3Cp. However, there 
is little else to guide one’s choice of integration constraints without anticipating the subsequent 
problem of relofting S to be a stream surface. 




Two Dimensions 


Figure 4, — Schematic of Tangential Vector Field on S 



In this report we do not address the theory and mechanics of stream surface lofting al- 
though specific examples will be described in section 6. For present purposes, we assume 
that every point of inflow on the perimeter of S remains fixed and every other point is dis- 
placed in the direction of the local normal (or in some other nonsingular direction). In 
general, S will not close, i.e., the points of outflow will not remain fixed. In order to obtain 
such closure, we replace the specification of at every point of inflow by the auxiliary 
condition 



3n 


de = o 


L 


( 11 ) 


where L denotes a streamline of t on S. The condition (11) attempts to force streamlines 
of the resultant flow passing through points of inflow to also pass through points of outflow. 


An alternative auxiliary condition can be employed when closure is unimportant, but it is 
desired that the resultant flow behave smoothly at each point of inflow to S (e.g. where S is 
a vortex sheet emanating from a wing leading or trailing edge). Here the available degrees 
of freedom can be used to annihilate the highest order singularity in the flow at such points. 
This can be achieved by specifying 


3 $ 

— = finite 
3n 


( 12 ) 


which in effect controls certain singularity strength continuity properties (see app. C). 

The use of the auxiliary conditions (1 1) and (12) together with equation (9) results in boun- 
dary values on S which are no longer exactly equivalent to the specification of potential, 
but from experience appear to be generally well posed. The major exception here is the case 
where D is finite (fig. 2) and potential is nowhere specified on B. 

We close this section by noting an important generalization of the boundary conditions of 
this section in the case where S bounds D on both sides (fig. 5). Then the specification of 
^ on one side may be replaced by the specification of A0 (the jump in <f> across S). Simi- 

b<f> 

larly, the specification of — on one side may be replaced by the specification of ^ — 

dn 

90 

(the jump in — across S). 

9n 
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4.2 DERIVATION OF INTEGRAL EQUATIONS 

Under rather general assumptions (ref. 1 , page 22 1 ), Green’s third identity shows that any 
solution of (1) may be expressed as the perturbation potential induced by a combination of 
source singularities of strength o and doublet singularities of strength pi distributed on the 
boundary B, i.e,, 


B B 

0 

Here r is the distance from the field point P to the boundary point Q and - — is the 

dnq 

derivative in the direction of the inner surface normal, i.e., directed into the domain D (see 
fig. 1 ). (On portions of B where both sides bound D the integration is performed over one 
side only.) If D is infinite, the representation (13) presupposes a certain behavior at infin- 
ity, namely that ^(P) is arbitrarily small when P is sufficiently distant from B. 
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Equation (13) represents an explicit solution to any boundary value problem of the previous 
section once densities a and fji have been determined for which 0 takes on the specified 
boundary values. For this purpose, equation (13) can be differentiated to yield 

B B 

(The underlined terms in equations (13) and (14) are generally referred to as kernel func- 
tions for the perturbation potential and velocity, respectively.) 

Upon sending P to B in (13) and (14) and substituting the right sides of these equations 
into the boundary condition equations of the previous section, we obtain integral equations 
for a and If the boundary value problem in D is well posed, these equations are suffi- 
cient to determine unique solutions a and fi on any portion B' of B where both sides 
bound D (fig. 6). To see this, we note (ref. 1, page 221) that 

/i(Q) = A0(Q) QeB' (15) 


and 


o(Q) = A 


90(Q) 

3n 


QeB' 


(16) 


where the symbol A denotes the difference between the limiting value of the quantity on 
the side of B' whose normal is n and the limiting value on the other side. However, on any 
portion B' of B which bounds an adjacent domain D' (fig. 7), an infinite number of dif- 
ferent source and doublet distributions can produce the desired potential in D. To under- 
stand this phenomenon, we note that equation (13) defines a potential in D' as well as D. 
Applying equations (15) and (16) in the present context, it follows that o and fx are unique- 
ly determined only when potentials in both D and D' are specified. Thus, the determination 
of source and doublet distributions solving the boundary value problem in D requires the 
specification of potential in every domain adjacent to D. This can be done explicitly, 
by specifying a known potential or implicitly by assigning boundary values determining a 
unique potential. For this purpose, the boundary B' of D' can be augmented by any num- 
ber of additional surfaces in the interior of D\ The choice of a potential here is guided by 
considerations of efficiency and accuracy. For example, it is advantageous to minimize the 
number of integral equations required to solve for a and fi by choosing the potential in D' 
so that either a or is specified on the boundary via equations (15) and (16). In this 
connection, the specification a = 0 or m = 0 reduces the computation cost in solving the 
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remaining integral equations. Finally, the potential in D' should be chosen so as to reduce 
numerical errors associated with large singularity strength gradients. 



Figure 6. — Portion of Boundary Which Bounds Domain on Both Sides 


A 



Figure 7. — Portion of Boundary Which Bounds Adjacent Domain 
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Several common choices illustrate these considerations. The first is the specification of zero 
perturbation potential in D\ From equations (15) and (16), we see that this choice implies 
that fjL and a are equal respectively, to the limiting values of 0 and ^ on the side of B' 

bounding D. Specification of ^ on this side is then equivalent to the direct specification 

dn ^ 

of o. This particular choice has the advantage that the potential in D does not by itself 
cause large gradients in a and fi, although it does not specifically reduce the gradients 
caused by the potential in D. The method of reference 1 1 has implemented a similar choice 
with considerable success. A typical boundary value problem for this choice is shown in figure 8. 


D 



Figure 8, — Zero Perturbation Potential in Adjacent Domain 


A second popular choice is the potential in D' with the same value on B' as the potential 
in D. This choice enjoys the advantage that g is zero on B\ since <f> is continuous across 
B'. The method of reference 2 pioneered this alternative. One problem with this choice in 
connection with lifting configurations is that the source strength becomes unbounded at any 
point of B' common to a surface B" bounding D across which potential is discontinuous, 
e.g.,awake. For this reason, the (doublet) surface B'' is usually continued into D' and 
assigned boundary values designed specifically to reduce the source strength gradients every- 
where on B'. (Depending on the specific technique, the solution of additional equations 
may be required.) The methods of reference 6 and 9 have used this device with excellent 
results. A typical boundary value problem for this choice is shown in figure 9. 
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A third successful alternative is the technique used in reference 8. Sources are again the 
primary singularity on B , but the requirement of additional boundaries in is avoided by 
allowing a limited (e.g., linear) variation of doublet strength on B'. A typical boundary value 
problem for this choice is shown in figure 10. 


bd) 



Figure 10. — Limited Variation of Doublet Strength on Boundary 
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A fourth alternative is the specification of a potential in D' with the same normal deriva- 
tive on B' as that of the potential in D. This choice enjoys the advantage that o is zero 
on B' since-^ is continuous across B'. Moreover, special treatment of intersecting doub- 
let surfaces is not required. However, if D' is finite, such a specification does not 

determine a unique potential in D' in view of our earlier discussion of the interior Neumann 
problem. This indeterminancy is reflected by the fact that a constant doublet density on a 
closed surface induces no velocity in the interior; therefore, it would be impossible to solve 
for a unique doublet density on B'. Here again, the problem can be remedied by introduc- 
ing an artificial boundary surface in D' on which potential is specified, e.g., a source sur- 
face on which </> = 0. If the specified values of M on B' satisfy equation (2), total source 

strength on such a surface will vanish. Otherwise, the source strength sum will be precisely 
that required to produce the net flow through B' desired. A typical boundary value prob- 
lem for this choice is shown in figure 1 1 . 



Figure 11. — Doublet Surface Modeling 



5.0 METHOD 


The present method uses a panel scheme to achieve a numerical solution to equations of 
section 4.2. Panel schemes proceed by first dividing the boundary surfaces into panels. 

In the present method, source (a) and/or doublet (/i) distributions are assigned to each 
panel. These distributions are expressed in terms of unknown singularity parameters Xj 
associated with the panel and neighboring panels. A finite set of control points (equal in 
number to the number of singularity parameters) is selected at which the boundary condi- 
tions are imposed. Evaluation of the boundary conditions results in a finite set of linear 
equations denoted symbolically by 


AX = b 


(17) 


Here, X is the vector of unknown singularity parameters, b the vector of specified boun- 
dary conditions, and A the matrix of “influence coefficients.’^ (Ay represents the influence 
of Xj on the boundary condition bj.) The solution of (17) for X may be accomplished by 
any number of computer algorithms; whereupon substitution of the corresponding distribu- 
tions of o and into equations (13) and (14) yields the potential and velocity at any 
point in D or on B. Pressure coefficients may be computed from the velocities; following 
which force and moment coefficients may be obtained by integration. 

The approach of the present method in performing these tasks is a building block approach 
in which the influence coefficient equation set (17) is defined by assembling independent 
networks of panels, each of which contributes as many equations as unknowns. A network 
is defined as a paneled portion of the boundary on which either a source or doublet distri- 
bution is defined accompanied by a properly posed set of analysis (Neumann) or design 
(Dirichlet) boundary conditions. All networks fall into four catagories: source/analysis, 
doublet/analysis, source/design, and doublet/design. The present method employs a variety 
of standard networks including one source/analysis, one doublet/analysis, two source/design, 
two doublet/design and two special wake/design networks. The construction of each network 
requires numerical development in three areas: A. Surface geometry definition; B. Singular- 
ity strength definition; and C. Control point selection and boundary condition specification. 
Essential features of the computational schemes in each of these areas are summarized below 
and discussed in detail in appendixes A, B, and C, respectively. 

A. Geometry input for a network is assumed to be a grid of comer point coordinates 
partitioning the network surface into panels. Panel surface shape is obtained by fitting 
a paraboloid to comer points in an immediate neighborhood by the method of least 
squares. (See app. A.) 

B. Discrete values of singularity strength at certain standard points on each network are 
assigned as singularity parameters (the Xj of equation (17)). Singularity splines are con- 
stmcted for each network type by fitting a linear (source) or quadratic (doublet) dis- 
tribution on each panel of the network to a subset of these discrete singularity param- 
eters by the method of least squares. (That is, the singularity distribution a(Q) or 
/i(Q) for each panel is expressed as a sum of products of singularity parameters in an 
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immediate neighborhood times linear or quadratic distributions of Q.) The least square 
technique allows minor discontinuities in value and gradient of the singularity strength 
across panel edges, although virtual continuity is assured if paneling is sufficiently dense. 
An alternate spline enforcing strict continuity of value at panel corners is available for 
the doublet/analysis network. For coarse paneling, this spline yields local doublet grad- 
dients which better reflect the global variations of doublet strength. (See app. B.) 

C. Certain standard points on each network are assigned as control points. These points 
include panel center points as well as network edge points. Boundary conditions in- 
volving the specification of potential or velocity are applied at panel center points for 
the purpose of controlling local properties of the flow, e.g., no flow through the sur- 
face. Auxiliary boundary conditions at edge control points serve to convey the proper 
continuity of singularity strength and gradient across network junctions, or else to con- 
trol global properties of the flow such as closure. In the case of design networks, these 
boundary conditions remove the degrees of freedom produced by specifying only tan- 
gential derivatives of the potential at panel center control points (see the discussion for 
equations (9), (10), (1 1), and (12)). In the process, they enforce standard aerodynamic 
auxiliary conditions such as the Kutta condition, closure condition, and the Helmholtz 
law. (See app. C.) 

Once each network has been constructed, the solution of the flow problem requires 
numerical development in three additional areas: D. Calculation of the influence co- 
efficients comprising the matrix A of equation (1 7); E. Solution of the matrix equa- 
tion (17) for the singularity parameters X; and F. Computation of aerodynamic quan- 
tities of interest. Essential features of the schemes used in each of these areas are sum- 
marized below and discussed in detail in appendixes D, E and F respectively. 

D. Two expansions of the induced potential and velocity kernels are employed. The near 
field expansion is based upon the assumption of a small panel curvature; the far field 
expansion requires a relatively large separation between the field point and panel. All 
resultant integrals are evaluated in closed form, using linear recursion relations that 
have as initial conditions the fundamental logarithm and arctangent transcendental terms 
appearing in flat panel, constant singularity strength techniques. (See app. D.) 

E. The Crout decomposition technique is employed by the pilot code to solve equation 
(17). This technique decomposes the matrix A into the product of a lower triangular 
matrix and an upper triangular matrix from which point the solution X is easily ob- 
tained by forward and backward substitution. Application of the technique is accom- 
plished with the aid of mass storage devices. For this purpose, the matrices (initially 
generated row by row) are partitioned into randomly accessible blocks. In-block partial 
pivoting is used to control error growth. (See app. E). 

F. Once the singularity parameters Xj are found and the corresponding source and doublet 
distributions determined, the potential and velocity at each panel center control point 
are calculated from equations (13) and (14). A distribution of pressure coefficients on 
each network is then found using the least square fit techniques of appendix B. 

Finally, this distribution is integrated to yield force and moment coefficients. (See 
app. F.) 
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6.0 RESULTS 


6.1 INTRODUCTION 

In this section, we present computed results demonstrating the numerical characteristics of 
the present method. These results were generated by a pilot computer program written in 
FORTRAN IV and COMPASS 3 languages for the CDC 6600 digital computer. Because of 
the versatility of the method, a complete analysis and display of all its capabilities would 
require an enormous amount of effort. Each application of the method, especially when 
design networks are involved, could well be the subject of a separate study and report, 
involving modeling alternatives, timing and convergence checks, experimental comparisons, 
etc. One such study has already been performed in connection with the application of a 
doublet/design network to the solution of leading-edge vortex flow problems (ref. 17). 

Hence, the results presented in this report are intended to give a broad view of the capabil- 
ities of the method. Because the method may be the basis of future production codes 
covering the subsonic, supersonic, steady and unsteady flow regimes, we have included cases 
which go somewhat beyond demonstrating mere feasibility. Such cases are the result of 
systematic efforts to test critical numerical features of the method, optimize modeling for 
accuracy and efficiency, and explore techniques for practical use. 

Numerical results are presented in sections 6.2 through 6.7 and are classified according to 
the network types featured. In section 6.2, we demonstrate use of the source/analysis net- 
work with some rather simple test cases. In section 6.3, we present results for the doublet/ 
analysis network. In contrast to source networks, the doublet networks have always re- 
quired more careful treatment, because doublets create a discontinuity in potential itself, 
and also induce a higher order singularity in the flow [equations (D.121) and (D.141)] . 
Therefore, we present cases testing convergence and sensitivity, related to the use of this 
network. Section 6.4 contains results for general analysis problems involving both source/ 
analysis and doublet/analysis networks. Cases are presented comparing alternate formula- 
tions. In section 6.5, we demonstrate the use of source/design networks with cases involving 
actual surface relofting. In section 6,6, we describe the application of a doublet/design net- 
work to the solution of leading-edge vortex flow problems (ref. 17). Finally, in section 6.7 
we discuss the numerical efficiency of the present method and pilot code. 

6,2 SOURCE/ ANALYSIS NETWORKS 
6,2.1 CIRCULAR CYLINDER 

An impermeable circular cylinder in a uniform flow was simulated with a type 1 (source/ana- 
lysis) network* containing 6, 12, and 18 equally spaced curved panels for the whole cylinder. 
The panels were elongated in the cross flow direction to approximate two dimensional flow. 


* Schematics of singularity parameter locations and control point locations are given in 
figures B. 1 . and C. 1 , respectively, for various network types. 
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Boundary conditions of type equation (C.3) with cy = 0, = 0 and Pjy~ 0 were applied 

at each panel center control point. Resultant source strengths and tangential velocity mag- 
nitudes are plotted in figure 1 2, along with exact data and results generated by earlier flat- 
panel, constant-strength methods. The 1 8 panel case is virtually indistinguishable from the 
exact solution. (In contrast, 35 flat, constant strength panels still produce a noticeable error 
in source strength. Due to fortuitous error cancellation (ref. 13), flat, constant strength 
source panels produce exact tangential velocities for a cylinder.) The six panel case begins to 
deviate significantly from the exact solution. Analysis shows the errors are primarily due to 
the violation of restriction (A. 2 1 ) of appendix A. The inadequate panel density results in an 
approximate surface which bulges to a radius of 1.15 at panel center points. 

6.2.2 SPHERE WITH RANDOM PANELING 

An impermeable sphere in a uniform flow was simulated with an 81 panel source/analysis 
network. For this purpose, advantage was taken of one plane of symmetry by paneling half 
the sphere only, and then calculating the perturbation potential as the sum of the potentials 
induced at a point and its image. A 10 x 10 comer point grid was generated, using a random 
number generator leading to a wide variation in panel size and shape as shown in figure 13a. 
The use of curved panels produced a geometry remarkably close to spherical. All 8 1 control 
points fell within a distance of 0.005 from the surface of the unit sphere, despite the fact that 
the maximum diameter of some panels spanned an arc of over 60°. As in the previous case, 
the use of such large panels violates restriction (A.21). (The consequences are not so serious 
in the present example because of the number and proximity of adjacent grid points used to 
obtain panel surface fits.) 

Velocity magnitude at each control point is plotted in figure 13b as a function of polar angle 
relative to the freestream direction. Agreement with the exact solution is good. In contrast, 
an analysis with flat constant strength source panels (representative of earlier technology) 
using the same panel arrangement, resulted in velocity magnitudes (V/Voo) that were scattered 
between 1.2 and 1.7 in the range 85 ° < 0 < 95°. This example demonstrates the extreme 
forgiveness of the method to irregular paneling, a feature which greatly enhances its practical 
usability for applications involving complex configurations where regular, evenly spaced 
paneling cannot always be arranged. 

6.3 DOUBLET/ ANALYSIS NETWORKS 
6.3.1 THIN CIRCULAR WING 

In this example, a thin circular wing at an angle of attack was simulated with a mean surface 
doublet/analysis network (type 2 with the least square doublet spline of section B.2 of appen- 
dix B). A doublet/wake network (type 8 with the least square spline) was abutted to the 
tr ailin g edge to yield a lifting solution. (See, for example, ref. 18, page 538 for the necessity 
of using a wake to generate a lifting solution.) Paneling for the right half of the wing and 
wake is displayed in figure 14a. On the wing, cosine spacing was employed along latitude and 
longitude lines with panels becoming triangular at the tip. The wake panels were attached to 
corresponding trailing edge panels and elongated in the stream direction. The entire wing and 
wake consisted of 108 effective panels (54 actual panels since symmetry was exploited). The 
freestream was directed in the plane of the wing and (linearized) boundary conditions of 
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(a) CIRCULAR WING AND WAKE PANELING 
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(b) SPANWISE LIFT DISTRIBUTION 
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Figure 14. — Circular Wing 
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type equation (C.3) with cy = 1 , Cl = 0, = - 1 employed to simulate unit angle of attack 

for comparison with angle of attack slope data. 

The resultant span wise circulation distribution is plotted in figure 14b. The plotted points 

were obtained from doublet strength along the trailing edge and agree well with the exact 

solution. Integration of this distribution produced a lift coefficient of 1 .776 versus the exact 

coefficient of 1 .790. (In particular, we have Cy = f See page 546 of Ref 18.) 

V^Sr J 

Pressures at three span stations are shown in figure 14c and also agree well with the 

exact solution. 

It is instructive to consider the role of the network edge control points in this case (see fig. C.l 
and section C.3 of app. C). The boundary conditions at leading edge control points automa- 
tically forced doublet strength to be zero on that edge since doublet strength was zero ahead 
of the leading edge. The boundary conditions at centerline control points forced the span- 
wise derivative of doublet strength to be zero on the centerline because of the implied presence 
of an image doublet surface. Finally, the boundary conditions at trailing edge control points 
in conjunction with those at opposing wake control points produced continuity of doublet 
strength and its transverse derivative onto the wake. Because of the particular construction of 
the doublet distribution on the wake network, this resulted in a vanishing chordwise deriva- 
tive of doublet strength at the wing trailing edge. Recalling that doublet strength is identical 
to potential jump (equation 1 5), we see that all the usual planar wing edge conditions includ- 
ing the Kutta condition were automatically satisfied. 

6.3.2 THIN SWEPT WING WITH RANDOM PANELING 

Figure 1 5 shows the stability of the mean surface doublet/analysis network (type 2 with least 
square doublet spline) under extreme conditions of panel size, shape and control point loca- 
tion. Here, 48 panels were used to represent the right half of a thin swept wing at 5.7° angle 
of attack. (The wake network is not shown but is similar to that of the previous case). Boun- 
dary conditions were specified as in the previous example. The wing panel layout was defined 
by means of a random number generator, resulting in panels that varied considerably in shape 
and size, that were occasionaly nonconvex, overlapping, and even inverted. Nevertheless, 
the spanwise lift distribution as shown in figure 15b is highly accurate. Chordwise pressure 
distributions displayed in figure 15c are likewise stable and deviate appreciably from the 
reference solution only near the leading edge where pressure becomes singular. Mismatches 
in doublet strength and derivative across panel edges occurred in this region, indicating that 
finer leading edge paneling is required for accuracy. (For the two values of tj shown, there 
are only two panels over approximately the first one-third of the wing chord at 17 = 0 and 
three panels at 1? = 0.5). 

6.3.3 THIN RECTANGULAR WING WITH VARYING PANEL DENSITIES 

For this case, a thin rectangular wing of aspect ratio 2 at 0.0 1 radians angle of attack was 
analyzed, using varying panel densities to check convergence of the solution. Symmetry 
was exploited and the right half of the wing was simulated with a mean surface doublet/ana- 
lysis network (type 2 with the continuous doublet spline of section B.3 of appendix B). The 
right half of the wake was simulated with a doublet/wake network (type 8 with the continu- 
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ous doublet spline). Both uniform and cosine spacing were employed with typical panel lay- 
outs shown in figure 16a. Cases were run with 4, 16, 49, 64, and 256 wing panels. Boun- 
dary conditions were imposed in the same manner as in the two previous examples. 

Resultant lift coefficients are plotted in figure 16b. (AU values have been scaled to one radian 
angle of attack.) Lift coefficients shown were obtained by integrating trailing edge doublet 
strength as before. Lift coefficients can also be calculated by integrating panel pressures as 
described in appendix F. For the continuous doublet/analysis spline, these values approach 
close agreement when paneling becomes sufficiently dense to produce negligible truncation 
errors. For example, the lift coefficients calculated by integrating panel pressures differed 
by 4.3%, 1.3%, 0.49%, 0.36%, -0.03% for the 4, 16, 49, 64 and 256 panel cases with cosine 
spacing, and 4.3%, 1.5%, 0.66%, 0.53%, -0.04%, for the same cases with uniform spacing. 

(Our observation has been that the deterioration of computed lift coefficients with decreas- 
ing panel density is less for values calculated by integrating trailing edge doublet strength, 
although in this case, the values calculated by integrating pressures happen to be better 
because of fortuitous truncation errors, due to rapid span variations near the tip.) The lift 
coefficients are apparently convergent with that of the 256 panel case using cosine spacing 
differing from that of the highly accurate and converged solution of reference 20 by less 
than 0.05%. The triangular symbols indicate the number of pressure modes rather than 
panels. The equal spacing cases converge somewhat slower than cosine pacing cases and 
this is probably a result of the relatively sparse paneling near the leading edge where pressure 
gradients are large. 

The locations of the centers of pressure are plotted in figure 16c and demonstrate conver- 
gence as well. Spanwise lift distributions are plotted in figure 16d and 16e. Chordwise 
pressure distributions at several span stations are plotted in figures 16f and 16g. These plots 
confirm convergence. The 16 panel case with cosine spacing is already sufficiently accurate 
for most purposes. The 64 panel case with cosine spacing is highly accurate. 

6.3.4 THIN RECTANGULAR WING WITH PANEL MISMATCHES 

The ability to refine panel density in regions of interest without refining panel density else- 
where is essential for the practical usability of any panel method. The present example 
demonstrates the characteristics of the present method in this regard, using the wing and 
network types of the previous example. In figure 17a, we display panelings for the right 
half of the wing and associated lift and moment data. In each case, three networks were 
used to represent the wing and two were used to represent the wake as shown on the left of 
the figure. Networks I, II and III were assigned panelings corresponding to the representation 
of the wing by 144 uniformly spaced panels. Network IV was assigned N x N, uniformly 
spaced panels with N = 2, 3, 4, 6, 12 and network V was paneled accordingly. 


It was anticipated that the only successful cases would be those with N = 2 and N = 6 where 
the edge control points of network IV were directly opposed by those of the adjacent net- 
works (see fig. C.l. of app. C.) Surprisingly, all cases displayed basic stability. The 
matching of pressure (ACp) and circulation (F) across the edges of network IV was virtuaUy 
exact for the case N = 6. Matching for the case N = 2 is displayed in figures 1 7b and 1 7c, and 
matching for the case N = 4 is displayed in figures 17 d and 17e. The circulation is of prime 
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(d) MISMATCHES ACROSS NETWORK LEADING EDGE (N =■ 4) 
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importance to the global properties of the flow and in all cases, it is virtually continuous 
across the network junctions. The pressure mismatch across the edge 0.5 is small in 

both cases. The pressure mismatch across the leading edge of network IV is small in the 

case of N = 2, although there is a slight anomaly at about 0.3. The pressure mis- 

b 

match across this edge in the case N = 4 shows the clear lack of opposing control points. 

Here, continuity of pressure is not enforced at any point on the edge and the only stabilizing 
influence is from adjacent panel center control points. The discontinuity appears to be con- 
fined to the network edge, but until more data is available, it is recommended that density 
refinements be made so that all edge control points of the coarsely paneled networks be 
opposed by edge control points adjacent finely paneled networks. 

6.3.5 THIN RECTANGULAR WING WITH VARYING SPANWISE GAPS 

In figure 18, we display spanwise lift distributions for a square wing separated from its image 
by varying spanwise gaps. A 10 x 10 uniformly spaced paneling arrangement was employed 
and symmetry was exploited. Tlie purpose here was to examine the ability of the present 
method (and in particular, the edge control point boundary conditions) to account for the 
nonuniform convergence of the lift distribution as the right and left halves of an .*^2 rectan- 
gular wing were brought together. 

At a separation distance greater than 10, each half was for all practical purposes, an isolated 
square wing as evidenced by the symmetric spanwise lift distribution. As the separation dis- 
tance decreased, the presence of the image wing was felt and the lift distribution away from 
the inboard edge began approaching that of an aspect ratio 2 wing. However, as long as the 
method perceived a gap, the load at the inboard edge continued to vanish in confirmity with 
the requirements of potential flow. At a separation distance of 10^, the method was unable 
to continue the limiting process, probably because the inboard paneling was simply too coarse 
to account for the severe local behavior. At a separation distance of 10"^®, the method believes 
the halves to be joined, but the edge control point boundary condition equations (section C.3) 
of appendix C) are still computed from the actual separated geometry resulting in some ano- 
malies. At a separation distance of 10"^^, these anomalies disappear and the lift distribution 
of figure 16d is achieved. This example illustrates that in potential flow, a gap which is neg- 
ligible physically, can have an enormous effect on the solution. The present method seems 
to account for such an effect, although numerical limitations exist for very small gaps. 

6.3.6 THINT-TAIL 

In this example, a thin T-Tail configuration at 0.01 radians angle of attack and yaw was ana- 
lyzed to test the functioning of the edge control point boundary conditions for nonplanar con- 
figurations. A square wing was used for each of the three components of the T-Tail. Each 
wing was simulated with a doublet/analysis network (type 2 with continuous spline) to which 
was abutted a doublet/wake network (type 8 with continuous spline). Since T-Tail compari- 
son data were unavailable, two different panel configurations were run, the first with 25 
uniformly spaced panels per component and the second with 81 uniformly spaced panels per 
component. Panel center boundary conditions were of type equation (C.3) with cy = 1, 

Cl = 0 and = 0. The resultant spanwise pressure distributions at half chord are displayed 
in figure 19a and the spanwise load distributions together with force and moment coeffi- 
cients are shown in figure 19b. Agreement between the data of the coarse and fine panelings 
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(b) TOTAL SPANWISE LOAD DISTRIBUTION 
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Figure 19. — (Concluded) 
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is good. The pressure and load distributions have the right character. In particular, the dis- 
continuity in ACp across the junction of the horizontal components is matched by the 
ACp of the vertical component at the junction. (This follows from the definition of ACp 
and the fact that the upper surface Cp distribution is continuous across the junction of the 
horizontal component.) The same property is true of the total load distribution. At present, 
leading edge suction terms are not included in force and moment calculation; hence, only 
those coefficients not significantly affected are displayed. The magnitudes and signs appear, 
in our judgment, to be correct. 


6.4 COMBINED SOURCE/ ANALYSIS AND DOUBLET/ANALYSIS NETWORKS 
6.4.1 WING-BODY ANALYSIS 

In figure 20, we present aerodynamic data for a symmetric wing-body configuration at 10® 
angle of attack. The fuselage is a body of revolution of thickness ratio 0. 1 1 . The wing is 
symmetric, 10% thick, and of an aspect ratio 5.6 with a leading edge sweep of 47®. The con- 
figuration was first analyzed by the method of reference 6, using 936 flat, constant-strength 
source panels on the standard wing and body surfaces, accompanied by 1 2 lifting elements. 
(This represents a typical number of panels used for wing/body applications with this method.) 

The paneling employed by the present method is depicted in figure 20a and comprises 160 
surface panels. The first such network contained all body panels forward of the wing. The 
second and third networks contained all body panels above and below the wing, respectively. 
The fourth network contained all body panels aft of the wing. The body was represented by 
four source/analysis networks with a total of 96 panels. The wing surface was represented 
by a 64 panel source/analysis network as shown. An internal lifting system (not shown) was 
used to create lift. The lifting system was composed of four networks. The first was a 32 
panel doublet/analysis network on the camber surface of the wing with stream surface (im- 
permeability) boundary conditions. The second was a 4 panel type 8 (doublet/wake. No. 1 ) 
network, emanating from the wing trailing edge. The third was an eight panel type 8 (doub- 
let/wake No. 1 ) network inside the body, extending the internal lifting system to the center- 
line. The fourth was a one panel type 10 (doublet/wake No. 2) network extending the 
trailing edge wake to the centerline. The continuous spline was used for all four lifting sys- 
tem networks. 

Spanwise load distribution comparisons are plotted in figure 20b and chordwise pressure dis- 
tribution comparisons are plotted in figure 20c. In addition lift, moment, and drag coeffi- 
cients are compared in figure 20b. All three coefficients agree well with the reference values; 
however, the drag coefficient agreement must be considered fortuitous in view of the extreme- 
ly sparse wing leading and trailing edge paneling. The calculation of accurate drag coefficients 
by integrating panel pressures generally requires a greater concentration of panels near the 
leading edge and trailing edge than the calculation of accurate lift and moment coefficients. 

It should also be pointed out that the drag values computed by the method of reference 6 
should not be considered as a valid standard, since that method has never been shown to pro- 
duce reliable drag values from integrated surface pressures. 
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(b) SPANWISE LOAD DISTRIBUTION 
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Figure 20. — Wing Body Analysis 




The spanwise load distribution agreement is good, although there is a slight discrepancy in 
average body lift. (The average body load was obtained by subtracting the total wing load 
from the total configuration load. The result was divided by the total load and distributed 
uniformly over the body span fraction. No precise span loading distribution for the body 
could be calculated because the body panel columns were not located along constant span 
stations.) Wing pressure agreement is excellent at the three inboard span stations. The dis- 
agreement in lower surface pressures at the outboard span stations is due to higher span- 
wise velocity components predicted by the present method. The reference method is known 
to underestimate spanwise velocities near wing tips, However, the discrepancy may also be 
due in part, to the width of the outboard panels employed in the present analysis. 

6,4.2 THICK WING ANALYSIS WITH FIVE DIFFERING 
BOUNDARY VALUE PROBLEM FORMULATIONS 

In section 4.2, we described a variety of boundary value formulations for solving a flow prob- 
lem where part of the boundary of the fluid domain also bounded an adjacent domain. It 
was noted there that an infinite number of different source and doublet distributions on 
this part of the boundary could all produce the desired flow in the fluid domain. In this 
example, we show results generated by the present method, illustrating this point. Speci- 
fically, we analyze the wing of figure 21(a) in five different ways. The wing has an iR 2 rec- 
tangular panform and a symmetric, 1 1% thick Boeing TR17 airfoil section. A 0. 1 radian angle 
of attack is assumed. Analysis of the wing by the method of reference 6 using 624 constant 
strength panels on the wing surface and 144 constant strength doublet panels for the internal 
lifting system produced the reference solution data of figures 21(b) and 21(c). The wing tip was 
left open for the reference calculation, an aspect that can affect in an unpredictable manner, 
the pressures near the tip, due to the possibility of inflow or outflow from the tip. The mo- 
ment coefficient Cjvjy is calculated about the leading edge. 

The panel layout used for the five applications of the present method is displayed in figure 
21(a). Each wing network employed in these applications comprises the lower, upper or cam- 
ber wing surface and contains four panel columns spanwise and 12 panel rows chordwise. 

The tip was paneled in only one instance and a network containing two panel columns and 1 2 
panel rows was used for this purpose. A four panel type 8 network was used to represent the 
wake in each case. The flve methods of analysis are described in the following five paragraphs. 

1 . The first analysis performed was similar to that of the reference method. Source/anal- 
ysis networks were placed on the wing upper and lower surfaces and a doublet/analysis 
network was placed on the wing camber surface (fig. 9). All three wing networks were 
assigned zero normal flow panel center boundar\^ conditions. For each source network, 
these conditions were applied on the side of the surface exposed to the exterior flow. 
Resultant aerodynamic data is compared to that of the reference method in figure 22. 

The data generally agrees well. As a check of the force coefficients of the reference 
method, this analysis was again performed with twice as many panels, resulting in values 

Cl = 0.261, Cmy = -0.0549, Cd 0.0128. 
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2. The second analysis performed used a formulation similar to that in reference 1 1. Both 

source and doublet/analysis networks were placed on the wing upper and lower sur- 
faces, (fig. 8), Boundary conditions requiring zero perturbation potential on the 
interior side of the wing surface were assigned to doublet panel control points (type 
equation (C.l 1) with ky = 0, kL = 1 and = 0). Boundary conditions requiring 
source strength to be equal to the negative of the normal component of freestream 
velocity were assigned to source panel control points (type equation (C.3) with c^ = 0, 
cj) = 1 , jjn “ The latter boundary conditions in conjunction with zero per- 

turbation flow in the wing interior cause the wing to be a stream surface in the exterior 
flow. Because of the crude leading edge paneling and to a lesser extent, the open tip, 
the specification of zero perturbation potential at doublet control points did not pro- 
duce precisely zero perturbation flow in the wing interior near the leading and tip 
edges. Consequently, the wing was not quite a stream surface in the exterior flow. 

This did not seem to make a great deal of difference, however, since most of the result- 
ant aerodynamic data (fig. 23) agrees well with the reference data. The disagreement 
in drag and moment coefficients may be due to a slight error in the way the pilot pro- 
gram calculates pressure coefficients when source and doublet panels are superimposed 
(see app. F). This error would affect primarily leading edge pressures. 

3. The third analysis performed reversed the formulation of the first analysis placing 
doublet/analysis networks on the upper and lower wing surfaces and a source analysis 
network on the camber surface (fig. 1 1). Zero normal flow panel center boundary con- 
ditions were applied to the doublet networks. Boundary conditions requiring the per- 
turbation potential to be zero were applied at all control points of the source/analysis 
network. Resultant aerodynamic data are displayed in figure 24. As an interesting note, 
the source strength on the internal source/analysis network turned out to be equal to 
the slope of the wing thickness. 

4. The fourth analysis performed used doublet/analysis networks only. Doublet/analysis 
networks were placed on the wing upper and lower surfaces and assigned zero normal 
flow panel center boundary conditions. The resultant spanwise circulation data were 
excellent, but the resultant pressure distributions were grossly in error and could only 
be improved by increasing panel density. Indirect boundary conditions similar to those 
of the second analysis were then substituted. Here, zero total potential on the interior 
side of the wing surface was specified at panel center control points in an attempt to 
force zero total velocity (stagnation) in the wing interior. (This would, in turn, cause 
the wing to be a stream surface in the external flow because normal velocity is contin- 
uous across a doublet surface as proven on page 1 70 of ref. 1 .) Stagnation was achieved 
near the centerline, but rather large velocities were present outboard and these were elim- 
inated by paneling the tip. Stagnant flow inside the wing implies that the total velocity 
on the exterior wing surface is equal to the surface gradient of the doublet strength (see 
eq. (C.IO)). Although stagnation was not quite achieved in the wing interior, the ex- 
terior velocities were calculated from doublet strength gradient anyway, resulting in 
slightly more, accurate values of Cp than those produced by actual exterior surface velo- 
cities calculated in the usual way. Resultant aerodynamic data are displayed in figure 25. 
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Figure 23. — /H 2 Thick Wing Analysis #2 
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Figure 24. - M 2 Thick Wing Analysis #3 
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Figure 25. — 2 Thick Wing Analysis #4 
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5. The fifth analysis performed used linearized boundary conditions and linearized geom- 
etry (chapter 17, ref. 1 8). Source/analysis and doublet/analysis networks were 
placed on the wing camber surface, which in this case was the x-y plane. As in the 
third analysis, boundary conditions requiring source strength to be equal to the slope 
of the wing thickness were applied to all control points of the source/analysis network. 
Boundary conditions of type equations (C.7) with c^^ = 1 and Cjj = 0 were assigned to 
each doublet panel control point. To be consistent with linear theory, the freestream 
was directed along the planform centerline and angle of attack effects were achieved 
by specifying /?„ = -0.1 at these control points. In addition, pressure coefficients were 

calculated using the linearized formulaCp = ’2 ^ • This resulted in values of ACp 
which were unaffected by the presence of the source network. Hence, lift and moment 
coefficients as well as span wise circulation (fig. 26) do not reflect thickness effects and 
agree with the data of figure 16. The primary effect of the source network in this case 
was to change the average of the upper and lower surface values of Cp to reflect thickness. 

For a pure thin wing, the upper and lower surface values of Cp would have equal magnitudes 
and opposite signs. Tliere is not a great deal of difference in the results produced by the 
first four analyses and without further study, it is difficult to favor any particular formulation. 
In fact, the choice of a formulation could very well depend on the specific application at hand. 
For example, a requirement for surface streamline tracing would clearly favor the use of 
formulations No. 2 and No. 4, where potential and velocities can be computed everywhere on 
the surface from the value and gradient of doublet strength. (Incidentally, this latter fact 
would allow somewhat more accurate integration of force and moment coefficients 
since the effect of the variation of pressure within a panel may be included in the cal- 
culation.) For extremely coarse panelings formulation. No. 1 is probably the most 
stable and accurate method of analysis, but it requires extra boundary conditions in- 
side a lifting body. (However, the internal lifting system can be quite crude. For most 
wings, the use of only three panel rows chordwise does not appreciably degrade results.) 

6.5 SOURCE/DESIGN NETWORKS 
6.5.1 DESIGN OF ARBITRARY AIRFOIL 

Figure 27 shows an application of the type 5 (source/design No. 2) network to a two- 
dimensional airfoil design problem. A NACA 65-010 symmetric airfoil at zero angle of at- 
tack was chosen as the nominal configuration. The arbitrary problem selected was a redesign 
of the airfoil between 20% and 90% chord producing zero Cp there. Analysis of the NACA 
65-010 airfoil was accomplished by means of three source/analysis networks placed between 
0% and 20% chord, 20% and 90% chord and 90%> and 100% chords, respectively, as shown in 
figure 27a. (Panels were elongated in the crossflow direction to approximate two-dimensional 
flow.) The resultant pressure distribution is displayed in figure 27b and is virtually identical 
to that given in reference 2 1 . The center network was then replaced by the type 5 (source/ 
design No. 2) network with tangential velocities of freestream magnitude as boundary condi- 
tions. Together with the closure condition (eq. 11) these boundary conditions produced a 
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flow with nonzero normal velocities at control points of the design network. The network 
panels were then relofted to eliminate the normal flow. This was done by sequentially moving 
each point a suitable distance in the positive or negative z direction, starting with the second 
point of the center network and continuing aft. The distance is determined by the require- 
ment that the unit normal of each updated panel be orthogonal to the velocity vector computed 
at the panel center prior to the update. In theory, the closure condition will ensure that the 
last point of the network will not be moved at all, thereby, maintaining continuity in geometry 
between the center design network and the aft analysis network. Because of the numerical 
discretization involved in the procedure, a slight gap may appear, in which case the whole 
design network is rotated about the initial comer point to achieve exact closure. Three itera- 
tions produced the reasonably converged geometry displayed in figure 27a, Analysis of this 
geometry produced the pressure distribution in figure 27b, which is close to that specified. 

6.5.2 REDESIGN OF TR 17 AIRFOIL 

A similar, but more practical application of the type 5 (source/design No. 2) network is 
shown in figure 28. Here, a Boeing TR17 symmetric airfoil (11% thick) was modified to pro- 
duce a different aft loading at 0° angle of attack. An analysis of the original airfoil section 
produced the pressure distribution plotted in figure 28b. The desired upper surface pressure 
distribution aft of 50% chord is also shown. In order to avoid the sharp discontinuities in 
pressure of the previous example, the upper surface geometry was allowed to vary aft of 30% 
chord and the original pressure distribution was specified between 30% and 50% chord. The 
design procedure used was identical to that of the previous example with the exception of 
the addition of an internal (doublet/analysis) lifting system and wake to allow for a lifting 
solution. Analysis of the designed section after one iteration indicated convergence. 

6.5.3 FULL AIRFOIL DESIGN 

Figure 29 shows a test of the type 5 (source/design No. 2) network on a full airfoil design 
problem. A circle was arbitrarily selected as a nominal airfoil. A 34 panel source/design 
No. 2 network was placed on the airfoil surface and assigned tangential velocities of an 18% 
thick, 3% cambered Karman-Trefftz airfoil as panel center boundary conditions. The closure 
condition (eq. 1 1) was also enforced and a lifting system composed of the usual internal 
doublet/analysis No. 1 network was employed to account for possible lift. The resultant flow 
produced nonzero normal velocities at source panel control points and these components 
were eliminated by relofting the airfoil panels in the same manner as the example in paragraph 
6.5.1 . The internal lifting system was relofted to the new camber line. The first iteration, 
shown in figure 29a, indicates the unsophisticated nature of the loft procedure used. No damp- 
ing was employed, i.e,, the geometry was allowed to assume the position predicted, even though 
the perturbation clearly violated linearity assumptions. The existence of a region of crossover at 
the leading edge did not cause the failure of subsequent iterations even though boundary con- 
ditions were applied on the interior side of the surface in this region. This is because sources 
alone were used to represent the airfoil surface and the specified tangential components of 
velocity automatically applied to both sides of the surface. At the fifth iteration, the cross- 
over was finally eliminated and the geometry rapidly converged. The airfoil produced on the 
ninth iteration is displayed in figure 29a and is almost indistinguishable from that of the true 
Karman-Trefftz airfoil. An analysis of this geometry with the source/design network replaced 
by a source/analysis network produced the Cp distribution shown in figure 29b which is 
close to that desired. 
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(a) GEOMETRY 
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(b) PRESSURE DISTRIBUTION 



Figure 29. — Full Airfoil Design 



6.5.4 WING DESIGN 


Figure 30 shows a three-dimensional test of a type 5 (source/design No. 2) network. Paneling 
of the wing-body model used for the test is displayed in figure 30a. For economy, the model 
is somewhat abbreviated ; nevertheless, the wing has camber, dihedral, and twist. The pur- 
pose of the test was to determine if the source/design network in conjunction with a geo- 
metry lofting routine could reproduce the original geometry from a perturbed geometry , 
using as boundary conditions the pressure distribution of the original wing. Analysis of the 
original geometry produced the wing pressures at four span stations displayed in figure 30b. 

The network arrangement and boundary conditions were identical to those of the example in 
section 6.4. 1 , except that the wing upper surface aft of the leading edge was simulated with a sep- 
arate source/analysis network. A perturbed wing geometry was achieved by compressing the 
wing upper surface approximately 30% at four span stations as shown in figure 30c. The con- 
figuration was reanalyzed, producing the modified wing pressures in figure 30b. The upper 
surface source/analysis network with modified geometry was then replaced with a type 5 
(source/design No. 2) network. The design boundary conditions ^pplied^to panel center 
control points of this network were of type equation (C.4) with tL = 0, ty = the unit vector 
in the direction of the modified geometry tangential velocities and Pi = the magnitudes of the 
original geometry velocities. A closure condition of type equation (C.17) with cy = 1, Cy = 0 
was applied to all four panel columns. The resultant flow produced velocities with nonzero 
normal components. The design sheet was then relofted to eliminate these components. For 
this purpose, all grid points at the design sheet edges except those at the wing tip were held 
fixed. The remaining grid points were allowed to move normal to the wing mean surfaces. 

The relofting was accomplished using the method of least squares, whereby a function was 
minimized with respect to a set of free parameters. The payoff was of the form 

PAYOFF = ^ 2 (V • n)2 

2 panel 
centers 


where ^ is the total velocity vector at a typical source/design network panel center, resulting 
from solution of the flow problem. The dependence of each panel center umt normal n on 
the four panel comer points is given in equation (A. 1 ). Let us assume th^ Pj is such^ cor- 
ner point, and that i is a unit vector normal to the wing n^an surface at Pj. "Oien if P: is one ; 
of the comer points allowed to move, we assume Pj(Xj) = Pj(0) + XjZ , where Pj(0) is the 
location of before relofting. Xj is then a parameter to be optimized in driving PAYOFF 
to a minimum. Many standard computational techniques are available for finding the mini- 
mum of PAYOFF with respect to all the parameters Xj. The particular technique employed 
in this case was a damped Newton— Raphson method (chapter 10 of ref. 22). 

Two iterations of this design procedure produced a geometry almost identical to the original 
geometry as shown in figure 30c. The difference between the first and second iterations was 
relatively small. Resultant pressure distributions for the designed geometry are displayed in 
figure 30b. At each lofting step, the payoff function was driven to a value very close to zero 
by unique values of the free parameters. This indicated that the closure conditions performed 
their function well. It also indicated the necessity of fixing the inboard edge grid points of the 
design network. 
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Figure 30, — Wing Design 





6.5.5 SPHERE DESIGN 


Figure 3 1 shows the application of a type 5 (source/design No. 2) network to an axisymmetric 
design problem. A sphere was analyzed with a type 1 (source/analysis) network. Symmetry 
was exploited so that only the right half of the sphere was modeled. Nine panel rows and 
columns were used for a total of 81 panels. The panel columns were defined by equally 
spaced lines of longitude with respect to poles at x = ±1 and the panel rows were defined by 
equally spaced x coordinates. 

Application of equation (C.3) with the usual impermeability option (cy = 1 , Cy ~ ~ 

produced tan gential ve locities agreeing well with exact velocities distribution 
V/Voo = 1 ,5.yi - X2 . The middle of the sphere was then chopped as shown in figure 3 1 a 
and the resultant geometry was analyzed with three source/analysis networks, the first com- 
prising the first two panel rows, the second comprising the middle five rows (chopped sec- 
tion), and the third comprising the last two panel rows. The resultant velocity distribution 
is displayed in figure 31b. The middle source/ analysis network was then replaced by a type 5 
source/design network with the tangential velocity distribution of the original sphere speci- 
fied at panel center control points. The closure condition (C.17) with cy = 1, cy = 0 was 
enforced on all panel columns. The resultant flow produced nonzero normal velocities at 
each panel center control point which were eliminated by the least square lofting technique 
of the previous example. In this case, grid points were allowed to move in a radial direction 
with respect to the x-axis. Only the front and aft edges of the design network were held 
fixed. Three iterations of this procedure produced a geometry indistinguishable from that of 
the original sphere. 
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(a) GEOMETRY <X-Y PLANE CUT) 
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Figure 31. — Sphere Design 
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6.6 DOUBLET/DESIGN NETWORKS 


6.6.1 LEADING EDGE SEPARATION 

Figure 32 shows the application of a type 4 (doublet/design No. 1) network to the analysis 
of separated flow off the leading edge of a sharp)-edged thin delta wing. This network was used 
to represent the free vortex sheet shed from the leading edge as shown in figure 32(a). The wing 
itself was modeled with a 30 panel type 1 (doublet/analysis) network. A simple kinematic 
model of the vortex core was simulated with a type 8 (doublet/wake No. 1) network. Other 
type 8 networks (not shown) were employed to represent wakes extending downstream 
from the wing and leading edge sheets. 

Analysis boundary conditions of type equation (C.3) with c^j = ^ ^n “ ® where speci- 

fied on the wing. A sophisticated iteration procedure was used to achieve the design-type 
condition of zero pressure jump across the sheet as well as the requirement that the sheet be 
a stream surface (refer to equations (5) and (8) of section 4.1). This procedure used a quasi- 
Newton scheme to simultaneously drive pressure and normal flow residuals to zero; however, 
at each iteration step, the effective design boundary conditions at free sheet panel centers 
were of type equation (C.8) with 





the subscript (0) denoting existing values. These conditions follow from the formula for Cp 
given by equations (F.l), frcm wl^h it is easy to show that the jump in Cp across a sheet ^ 

(ACp) is given by ACp = -2 • Vp)/V<>o^ where and Vp are given by equations (C.5) 

and (C.6) respectively. From equation (C. 1 0), we see that fw a doublet she^ 

ACp 2 . ■yj[x)/Voo^. Then -Vi/Voo^ (^Cp + ACp^) = (V^q • VMq) “ ^Ao> ~ VMo)> 

and the result follows by setting ACp = 0 and neglecting the last term on the right, because it 
is of second order in the changes from the existing values. Boundary conditions of type equa- 
tion (2) were specified at the control points along the free sheet edge adjoining the wing to 
enforce a Kutta condition there. 

Free sheet comer points (not attached to the wing) were allowed to vary in planes perpen- 
dicular to the centerline to satisfy the stream surface requirement. The variation was defined 
by fixing panel edge lengths in those planes and assigning slope angles as free parameters. 

Figure 32(b) shows detailed pressure distributions for a thin delta wing of aspect ratio 1.4559 
at a 14° angle of attack. Upper and lower surface pressure distributions agree with the experi- 
mental data from ref. 23 in spite of the sparse wing paneling. Figure 33 displays results for 
a thin delta wing of aspect ratio 1 . Normal force coefficients at four angles of attack are 
plotted in figure 33(a) and they agree well with the experimental data of Peckham (ref. 24) 
and theoretical results from the leading edge suction analogy of Polhamus (ref. 25). A 
typical load distribution is plotted in figure 33(b) and it agrees well with experimental data. 
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(a) PANELING OF DELTA WING AND VORTEX SHEET 



(b) SURFACE PRESSURE DISTRIBUTION OF DELTA WING 



Figure 32. — /R 1.4559 Delta Wing with Leading Edge Separation 
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(a) NORMAL FORCE COEFFICIENT AT VARIOUS a 



(b) TYPICAL LOAD DISTRIBUTION 



Figure 33. — 1 Delta Wing with Leading Edge Separation 
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6.7 NUMERICAL EFFICIENCY 


Enough information is currently available to surmise that the computing costs associated 
with the present use of curved panels and locally varying singularity strengths are quite com- 
parable on a panel by panel basis with those of first order panel methods. (In terms of 
total computing cost per problem the present method has clear advantages because of the 
reduced number of panels required.) 

Computing costs depend primarily on central processing (CP) time and mass storage use (10). 
Although 10 costs are important, they are little affected by the use of higher order panels 
and singularities. To be sure, additional geometry and singularity defining quantities must 
be stored for each panel, but 10 costs usually depend more on mass storage access requests 
than the quantity of data stored per record. For cases with less than a certain number of 
panels (about 750 in our case) the 10 costs of a well-optimized program are usually small 
compared to CP costs. For cases with greater than this number of panels, the 10 costs of 
solving the influence coefficient equations are dominant. 

For an influence coefficient program CP costs behave according to the rough formula, 

Cost = Cjn + C2n^ + C^n^ (18) 


where n is the number of panels. The linear term on the right is associated with the com- 
putation of panel geometry, singularity strength and control point defining quantities as 
well as the calculation of output. The quadratic term corresponds to the computation of 
the influence coefficients and the cubic term is associated with solving the influence coeffi- 
cient equation set. The severity of the cubic term can sometimes be reduced through the 
use of iterative techniques, but the diagonal dominance required for convergence is difficult 
to ensure for truly arbitrary configurations with interacting components. The quadratic 
term is somewhat too severe to represent the behavior of the CP time for generating influence 
coefficients when n is in the usual range of interest. This is because pure increases in panel 
density are accompanied by a higher proportion of efficient far field calculations. Neverthe- 
less, the quadratic assumption is sufficient for the present analysis. 

The coefficient Cj is considerably greater for the present method than for first order tech- 
niques. With the exception of output calculations, the linear terms dominate only when 
n < 25 and by the time n = 60, they contribute less than 10% of the total cost. However, 
for multiple freestream directions, the output calculations can represent up to 25% of the 
total CP costs even when n = 500. This is because the present pilot code calculates an enor- 
mous quantity of output— primarily for diagnostic purposes. Presumably, a production pro- 
gram would not require such extravagance. 

The coefficient C3 is unaffected by the use of higher order panels and singularities. Current 
projections indicate the cubic terms become dominant when n > 1000, so that the cost advantage 
of using a higher order method is clear for large cases. Here, however, we must note that the 
effective value of n increases for networks with edge control points. In most cases to which 
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current first order techniques are applied, the edge control point equations are somewhat 
uncoupled from the total equation set, although advantage is not derived from this fact in the 
pilot program. The reason for this uncoupling is that the edge control point boundary condi- 
tions effectively control singularity matching across network edges (app. c) and therefore, 
each such boundary condition equation depends only on a few unknown singularity para- 
meters in the neighborhood of the network edge. 

For problems with moderate numbers of panels (50 < n < 1000), the quadratic term of 
equation (18) is dominant. The calculation of influence coefficients is the area where some 
higher order methods would ordinarily be expected to incur substantial penalty. For the 
present method, this is not the case; primarily because the contributions of higher order terms 
can be expressed as simple combinations of lower order terms. Table 1 shows near field in- 
fluence coefficient calculation time comparisons between the present method and two first 
order methods— those of references 6 and 26. Table 2 shows analogous data for intermediate 
and far field calculations. Reformulation of the pilot code for the curved panel option has not 
yet been accomplished. The near field flat panel source calculation of the present method is 
faster than that of TEA 230, although computation of higher order terms requires about 10 
more operations. Actually CP time is only roughly proportional to operation count and such 
a variation falls well within the range expected for different codes. The same is true of the 
source intermediate and far field calculations. The times displayed in table 1 and 2 corres- 
pond to the computation of the potential and velocity at a field point induced by a singularity 
distribution on a panel. There is also an additional cost for transforming the resultant veloci- 
ties to global coordinates. Moreover, for the present method, there is a further cost of dis- 
tributing the influence coefficients to the singularity parameters. The net additional cost 
for the present method is approximately 0.24 milliseconds per source influence coefficient 
and 0.58 milliseconds per doublet influence coefficient. These costs dimmish the effective- 
ness of the far field computations in the present pilot code; however, they can be recovered 
for the important source/analysis network by implementing the third recommendation of 
section 7.2. 

In figure 34, we show CP time comparisons between the highly optimized TEA 230 program 
and the pilot code of the present method. The range of CP times for a given number of singu- 
larity parameters reflects variations due to source and doublet and to near field and far field 
ratios. TEA 230 appears to be about 30% faster on the average on a panel for panel basis. 

It is estimated that for a well-polished version of the present method, this advantage could be 
reduced to 15%. 
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Table /. — Near Field AiC Calculation Time Comparisons 


AIC 

Panel 

geometry 

CDC 6600 Time 
(milliseconds) 



Source 

Doublet 

1. Present method 

Flat 

8. 

15. 

(original formulation) 

Curved 

13. 

29. 

2. Present method 

Flat 

4.2 

6.3 

(reformulation of recursions, etc.) 

Curved 

6.1 

9. 

3. Present method 

Flat 

.91 

1.41 

(reformulation of code)^ 

Curved 

- 


4, TEA 230 (Ref. 6) 

Flat 

.94 

1., 14.* 

5. Flexstab (Ref.26) 

Flat 

Not 

easily 

isolated 

5.4** 


^Funded by Boeing IR & D 
* Linearly varying vortex 
** Constant strength vortex 


Table 2. — Intermediate / Far Field AIC Calculation Time Comparisons 


AIC 

Panel 

geometry 

CDC 6600 Time (milliseconds) 

Source 

Doublet 1 

Far 

field 

Int. 

field 

Far 

field 

Int. 

field 

1 . Present method 


.6 

.9 

.7 

1.8 

(original formulation) 


.6 

.9 

.7 

1.8 

2. Present method 


.09 

,54 

.26 

.72 

(reformulation of 


.09 

.54 

.26 

.72 

recursions, etc.) 






3. Present method 

Flat 

.09 

.54 

.26 

.72 

(reformulation of code)^ I 

Curved 

.09 

.54 

.26 

.72 

4. TEA 230 (Ref. 6) 

Flat 

.16 

.26 

Not comparable 

1 

5. Flexstab (Ref.26) 

Flat 


Not corr 

- - - -1 

iparable 



^Funded by Boeing IR & D 
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7.0 CONCLUDING REMARKS 


7.1 CONCLUSIONS 

An advanced panel method for solving boundary value problems associated with the analysis 
and design of arbitrary configurations in subsonic potential flow has been presented. Theory 
and computed results indicate the method has the necessary characteristics for widespread 
acceptance by the user community. These characteristics include the following: 

1 . Generality: The method is capable of solving a wide range of analysis and design boundary 
value problems. In particular, the method can handle virtually all analysis problems to 
which current panel techniques are applied. The added feature of combined analysis 

and design provides a capability that is extremely powerful for an extensive variety of 
applications. Most design problems involve the aerodynamic design of one or more 
components of a configuration in the presence of others whose geometrical shapes are 
fixed. The present method provides this capability along with the limiting cases of 
pure design or analysis. 

2. Flexibility: The method offers the user a variety of modeling options as well as a 
straightforward means of implementing those most suited to a specific application. For 
economical analysis and design, all the usual thin surface approximations are available. 

For more accurate results, exact boundary conditions may be applied on actual config- 
uration surfaces. Here, a variety of surface singularity distributions are available. The 
efficient source-alone option may be employed on nonlifting portions of the configura- 
tion. Components which shed vorticity can be modeled with combined source and doub- 
let surfaces or source surfaces accompanied by interior doublets or doublet surfaces 
accompanied by interior sources. The shed vorticity can be fixed in location and direc- 
tion or designed to satisfy pressure jump and stream surface requirements. These and 
other options are easily implemented through the use of standard networks to simulate 
arbitrary boundary surfaces. Aside from modeling versatility, the method offers sub- 
stantial paneling flexibility. No restrictions appear necessary on panel shape or orienta- 
tion so long as density is sufficient for the accuracy desired. 

3. Stability and Accuracy: For well-posed boundary value problems, the present method 
has been demonstrated to be numerically stable and accurate even under extremely adverse 
conditions. It displays a marked insensitivity to paneling layout and achieves acceptable 
accuracy with relatively sparse panel densities. Convergence to highly accurate results 
occurs at moderate panel densities. 

4. Efficiency: On a panel-by-panel basis, the present method appears capable of the 
same efficiency as that enjoyed by existing first order techniques. On a case-by-case 
basis, the present method has significant advantages because of the reduced number of 
panels required. These advantages become overwhelming for cases involving complex 
configurations. 
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7.2 RECOMMENDATIONS 


During the development of the present method, a number of improvements were discovered 
but not implemented because they did not sufficiently bear on establishing feasibility. None- 
theless, these improvements could have significant impact on any production-oriented ver- 
sion of the method and therefore we sketch the more important ones below. 

1 . Curvature expansion: The curvature expansion of equation (D.l 8) is valid for any field 
point, although it is used only for near field calculations. More efficient expansions are 
available for this case. As an example, we have 

1 1 

^ y^l^ - ha) J - (1 + ha)x ] ^ + [(1 - hb)T? - (1 + hb)y] “ + h“ 

when (x,y) is in close proximity to S. By a suitable linear transformation, this approxi- 
mation reduces to the flat panel expression with corresponding simplifications to equa- 
tions (D.22), (D.28), (D.32) and (D.34). Here, curvature is simulated by modifications 
to the panel shape and field point location. 

2. Near Field Influence Coefficient Calculation: One of the mitigating factors in the use of 
higher order singularity splines is that the complexity of the near field influence coeffi- 
cient calculation need not grow at the rate one might expect. This is because the in- 
creased continuity of the splines across panel edges allows the elimination of many can- 
celing line integrals. (Refer to sec. D.5 of app. D). At present, no advantage is derived 
from the fact. 

3. Projection Algorithm: For the construction of the higher order panel surface represen- 
tations and singularity distributions, an orthogonal tangent plane projection algorithm 
is currently used. (Refer to eq. (A.9) and (B.5), (B.2) and (B.3.) Because of the ten- 
dency of the user community to violate hypothesis (A.23), a length and azimuth pre- 
serving projection would be better. For example, the standard projection could be altered 
by scaling the distance from the projected point to the origin so that it is the same as 

the distance from the original point to the origin. 

4. Curvature: At the present time, violation of the curvature restriction (A.21 ) is ignored. 

It would be better to incorporate this constraint into the least square panel surface fit 
algorithm so that it is satisfied automatically. At the same time, provision should be 
made for better curvature definition in cases with sparse paneling, e.g., extra grid points 
at panel centers. 

5. Specification of Singularity Strength: The specification of source or doublet strength at 
a control point now requires an influence coefficient equation, although the equation 
usually amounts to direct specification of an element of the solution vector. Such an 
equation should be eliminated from the full equation set before solution. In the case 

of the continuous doublet/analysis spline wherein the values of the singularity param- 
eters do not directly correspond to singularity strength, the network boundary condi- 
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tion equations for any such network still uncouple from the full set and can be solved 
prior to the full solution. This reduction can also be applied to certain situations 
where doublet strength gradient is specified at panel center control points of a design 
type network. 

6. Superimposed Doublet and Source Panels: In the present numerics, superimposed doub- 
let and source networks are treated independently even when paneling is identical. 

There are clear advantages to identifying superimposed panels in the geometry defini- 
tion, influence coefficient calculation and force computation procedures. These advan- 
tages are best exploited by defining new network types with combined source and doub- 
let splines and duplicate control points. 

7. Force and Moment Calculations: The force and moment calculations of appendix F 
should be generalized to include momentum transfer terms for cases involving permeable 
surfaces and edge delta function terms for cases involving thin surface approximations. 

These generalizations are already under study in connection with the work of reference 17. 

8. Panel Center Boundary Condition Generalization: In connection with certain design 
problems, it appears desirable to require that the pressure at a given control point be 
equal to that at another control point. More generally, it appears desirable to allow 
boundary conditions which relate quantities at any number of control points. The 
closure condition (eq. (C.17) is already a boundary condition of this type. Another ex- 
ample is the slotted wind tunnel wall boundary condition, i.e., k^rf) + ^ = k^<l>n + 

^ 3n ^ T — 

3n 

where ~ evaluated upstream of the slotted wall section. The well-posed nature 

3n 

of the resultant boundary value problems is open to question, but preliminary experiments 
indicate that such boundary conditions can work very well in practice. The necessary 
modifications to equations (C.3) and (C.4) are straightforward. 

9. Grid Layout Generalization. The requirement that grid points and panels be arranged 
in rows and columns limits the ability to vary panel density in complex configurations. 
Clearly our least square definitions of panel geometry and singularity strength do not 
require such an arrangement and a more general panel identification logic is highly desirable. 

10. Design Procedure Generalization: Of all the design procedures described in the exam- 
ples of section 6, that of reference 17 appears to offer the most potential for practical 
design. Here, the problem of designing a stream surface to support a given pressure 
distribution is attacked as an optimization problem in which a payoff consisting of a 
weighted sum of residuals is minimized using sophisticated iteration techniques. The 
residuals include deviations in both pressure and normal flow and the optimizing variables 
are the singularity parameters as well as the parameters controlling surface perturbations. 
The pressure specification and stream surface lofting problems are coupled in this for- 
mulation, but this has the advantage of allowing trade studies. For example, inequality 
or equality constraints can be applied to the surface perturbation parameters, resulting 

in either an acceptable or unacceptable degradation in the specified pressure distribu- 
tion. Moreover, residuals of certain specified functions of pressure (forces and moments) 
can be weighted and added to the payoff function. 
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APPENDIX A 
GEOMETRY DEFINITION 


A,1 INTRODUCTION 

In this appendix, we shall define how a surface is constructed which approximates the true 
surface (i.e., analytically defined input surface) for the purpose of achieving a numerical 
solution to the flow problem. We assume that the overall configuration has been parti- 
tioned into networks. In the interior of each network, the true surface is assumed to have 
continuous position, slope and curvature. Any discontinuities in these quantities, as far as 
the true configuration surface is concerned, must therefore, occur at network edges. We 
assume that a discrete representation of the true network surface is provided by a grid ^f 
corner points r (I,J); I = 1,M and J = 1,N where the elements of the position vector of P are 
resolved in a global (x,y,z) coordinate system. A planar schematic of these points is shown 
in figure A. 1 , The grid layout is shown as rectangular for convenience only. For example, 
the wing comer points of figure 1 5a form an M = 9 by N = 7 grid. For each point J iden- 
tifies the column of points to which it belongs and I identifies the row. The grid is assumed 
to be sufficiently dense that the portion of the true network surface lying between four ad- 
jacent comer points does not deviate significantly from their average plane (to be discussed 
later). Moreover, we assume that the projection of the four corner points onto their average 
plane defines a quadrilaterial which at most degenerates into a triangle. (We, therefore, 
require at least three of every four adjacent corner points of a grid to be distinct.) The con- 
struction of an approximate surface or panel S lying between four adjacent corner points 
of the grid is the subject of the next two sections. 
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Figure A, 1. — Schematic of Corner Point Grid 
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A.2 FLAT PANEL APPROXIMATION 

Consider fouradja<^nt corner points P(I,J), P(I,J+1), P(I+1,J+1) and P(I+1,J) and relabel 
these points as?l, P2, ?3 and ?4 respectively. In this section, we wHl constnjct a first order 
approximation S to the true surface of the network lying between PI , P2, P3 and P4. At 
the same time, we will define a local orthogonal coordinate system on S for later use. 

A A ^ 

Let the orthogonal unit vectors 17 and f be defined as: 

. [P|tP4-P2-P3] 


. , [p] + P2 - P3 - P4J 
IPJ+P2-P3-P4JI 


A A A 

T?= f® f 


(A.l) 


where ® denotes the vector cross product. 
Let the average point Pq be defined by 


^0 = 4 


p, +P 2 + P 3 + P 4 


(A.2) 


The plane passing through Pq with normal ? is defined being the average plane of the points 
P^j, P9, ?3 and ?4. This plane is average in the sense that it also passes through the midpoints 
of the line segments joining ?i and ^2> ^2 ^3> ^3 ^nd^4, and P i as shown in figure A.2. 

To see this, we i^ote that it is sufficient to prove that f . (Pm -Po> = 0 where Pm is one of the 
midpoints, i.e.. Pm ~ (^i'i’^i+l)- direct substitition, one can show that for any i, 

“ ^o “ -^^1 + ^4 ” ^2 " ^3) 1 ^2 “ ^3 " ^4)- result then follows from (A.l), 

(A.2) and the well known fact that for any vector ? and b, "a . (?® b) = 0. 
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Figure A.2. - Average Plane of Four Adjacent Corner Points 

Note that the average plane of four arbitrary points is not unique and depends on the order 
assigned the points. Also note the cyclical order assigned the points above is such that the 
normals f of all the panels lie on the same side of the network surface. 

T[he points formed by the intersection of the average plane with the perpendiculars from P j , 
P 2 , P 3 and P 4 onto the average plane define a plane quadrilateral 2 whose boundary con- 
sists of the straight line segments joining these projections in cyclical order (see fig. A.2.). The 
panel 2 is taken to be the surface S which approximates the true surface. The coordinate 
directions of the local coordinate system associated with S ate taken to be the ^ and f 
directions, however, the origin is chosen as the centroid of the panel 2 , rather than the 
average point Pq. This choice makes little difference numerically, but is made to facilitate 
comparisons with earlier codes such as TEA 230, which also have origins at the panel centroid. 
The centroid P^ may be computed as follows: 

Define 


P5 = Pl 


(A3) 
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and 


= (Pi-Po)-i 

Tji = (Pj-Po).n 




Then 


4 

Pc = Po + ^ S Di [(^i + Ji+l) t+(r?i + T?i+i)^] 
i = 1 


where 


Di= yCSi T?i+i -TJi ^i+i) 


and 


4 

D. 2 Di 

i= 1 


(A.4) 


With this choice of local coordinates, any point P with components (x,y,z) in the global 
system is the same as the point Q with components (^, n, z) resolved in the local system 
if and only if 

Q = a(p-Pc) (A.5) 


i.e.. 


/i\ Mil Ai2 Ai3Yx-Xc\ 

I V J = l A21 A22 A23I y-yc I 

\f/ \A31 A32A33 /\z-Zc/ 

where (X(,, yc, Zc) global coordinates of the centroid P^ . 
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Here, A is the 3x3 orthogonal matrix, 


— \ 


(where , {r}] and are the column vectors containing the x, y and z components of 
the unit vectors and respectively, and T denotes transpose.) We also note the inverse 
transformation, i.e., 


P = aT Q + Pc (A.7) 

In general, the flat panel approximation to the network surface described above possesses 
discontinuities in position and slope along panel edges. However, because of our assumptions 
about the smoothness of the true network surface, these discontinuities become small rela- 
tive to panel size as the grid density increases. 

For a somewhat coarse grid, these discontinuities as well as deviations from the true surface 
can be large, relative to panel size; moreover, the effect of local curvature on the local flow is 
lost. For these reasons, a curved panel approximation has been devised and is described in 
the following section. 


A.3 CURVED PANEL APPROXIMATION 

In this section, we will construct a second order approximation to the true surface of the net- 
work lying between >*^2’ ^3 "^4* Our smoothness assumptions about the true network 

surface imply that for a sufficiently fine grid, the portion of the true surface lying between 
, ^2> ^3 ^4 approximated by a parabolic panel S whose shape depends only 

on’local grid points. The panel S may be defined in a variety of ways, but we select the 
following. Using the local (f, i?, f) coordinate system constructed in section (A.2), we assume 
S may be represented in the approximate form 

, ^) = f o + ^ ^ (A.8) 

The coefficients ( ro> > h' hb hv' frjr?) are^then obtained by fitting S to a set of points 
n which is composed of the corner points Pj P 3 and P 4 , as well as every grid point ad- 
jacent to these points. The fit is accomplished by the method of weighted least squares. To 
be specific, we minimize the quantity 
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where (^|^, t 7 j^, J'l^)are the known coordinates of a point in II. The weights Wj^ are 
chosen to be unity, unless happens to be a comer point of S, in which case Wj^ is 
chosen large (nominally 10^). With this choice, S very nearly passes through its comer 


poin^. If we let C be the column vector with components ( fr?’ fr/r? ) 

and Vj^he column vector with components ( 1, ^Vk^) , then f (£ki%) 

= Vi,T r, - ■ • - - - 


Vk^ C]^, and minimizing R with respect to the components of C yields 

(6x6) " “ (6x1) 


C = 
(6x 1) 


[s WkVkVjl 


I "kfkVk 

k 


(A. 10) 


(Here the superscript -1 denotes matrix inverse.) 

Once the coefficients in C have been obtained the expression (A.8) for the surface S can be 
simplified by a suitable transformation of coordinates. If the coefficient f is nonzero, the 
term may be eliminated by a rotation of the ({, tj, f) coordinate system about the f axis 
through the angle \{/e [-tt/4, tt/ 4] where 



J 

- 



— 


’ ' ^vv) ^ ® 

1 otherwise 


This defines a primed local coordinate system such that a point P having coordinates (x,y,z) 
in the global system is the same as a point ^ having coordinates (^', tj',?') in the primed 
local system if and only if 


Q = A'(p-Pc) (A.12) 

where 

a' =da 
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Here A and are defined in section (A.2) and 

^ cos\jj sini// 

D = ( - sini// COS)// 

,0 0 


In the primed local coordinate system S has the representation 

r=fo + dr+eT? '+3^2 + 57, -2 


(A.13) 


where 

a = cos2|// + cos\J(sin\J/+ sin2\//) 

b= sin2 1 ^- cosi/«ini/' cos^i/^) 

d = cos»// + fjp sini// 

e = sin0 + cos\jj 

If a =?^ 0, the term linear in ^ ' may be eliminated by a translation of the origin of the 
(^^ '>7^ coordinate system to the point 



r?' = 0, r = 0 


(A. 14) 


If b =5^ 0, the term linear in r?' may be eliminated in a similar fashion. If a = 0, then the term 
linear in may be eliminated by a rotation of the coordinate system about the r\ 

axis through an angle ae [-tt/2, tt/2] , where 



(A.15) 


If b = 0, the term linear in r} may be eliminated similarly by a rotation about the f axis. 

If a is nonzero but much smaller than d, equation (A. 14) translates the origin of the (f', 

T}\ f') coordinate system far away from the panel causing eventual numerical difficulties. 
The transformation (A.15) however, is a mere coordinate system rotation and clearly to be 
preferred. If this transformation is used when a 0, the expression for S in the resultant 
coordinate system will have no linear f'' term as desired, but will have quadratic 
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terms involving These quadratic terms are actually of higher order in 6 (see eq. (A. 21)) 
and can be ignored if the grid is sufficiently fine or removed iteratively otherwise. To re- 
move the quadratic terms involving iteratively we use the following procedure. We define 
a double primed local coordinate system in a manner similar to that suggested above, such 
that a point Shaving components (x,y,z) in the global system is the same as a point ^ 
having components 77", f") in the double primed local system if and only if 



K 


X-Xc 

Q = A"(P-Pc) or 


= [A"l 

y-Vc 


Irj 


Z-Zc 


(A. 16) 


where 


and 


and 


A = EA 


- (' ” 

E = I 0 cos ^ 

\ 0 - sin P 



0 

1 

0 


sin a 
0 

cos a 


/ cos ot 0 

I - sin j3 sin o; cos ^ 

\ - cos /3 sin a - sin j3 



a 

/3 


sin ^ ; 0 € [-tt/ 2 , tt/ 2] 

\\/i +d2 y 

sin"Y / ^ ; P € 1-7 t/ 2, 7 t/ 2] 

Vv/l +d2 + e2y 


(The matrix E represents the composition of the coordinate transformations for eliminating 
d and then the resultant e.) We then replace the original (f, 77, f) local coordinate system 
with the (f ", 17", f") system and repeat the construction from equation (A.8) on. If the pro- 
cedure fails to converge the flat panel approximation of section (A. 2) must be used. (This 
contingency has never been necessary.) For a sufficiently fine grid, the procedure will con- 
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verge rapidly so that shortly a coordinate system will be obtained in which d and e type terms 
in equation (A. 13) will be essentially zero and no quadratic terms in f will be present. The 
local coordinate system for which this occurs is then defined by the transformation of equa- 
tion (A.l 2) where now 


^ ( ^n^n-1 ^n-1 • ^1^1 ) ^ 


(A.17) 


and Dj, Ej are the D and E matrices for the ith iteration. 

The coefficient of equation (A. 13) may be eliminated by a translation for the f') 

coordinate system to the point 

= + (A.18) 


We then obtain a final local (^, 77 , f) coordinate system defined by the transformation from 
global to local coordinates, 


Q = a(p-Rq) (A.19) 

where A is now the matrix A', defined in (A. 17) and is defined by (A.18). In this 
coordinate system, S has the representation 

77 ) = + b77“ (A.20) 


The domain of ({, 77 ) here is again the quadralateral 2 defined by projecting the comer 
points ^3 and ^4 onto the (f , 7 ?) plane (see fig. D. 1 ). Equation (A.20) represents 

a considerable simplification over equation (A. 8 ) and leads to substantial economics in the 
calculation of influence coefficients (app. D). 

This completes the construction of the second order approximation to S. In section (A. l), 
we made the assumption that the grid was sufficiently dense that the true surface lying 
between four adjacent corner points did not deviate significantly from their average plane. 
This assumption implies that the coefficients a and b are limited in magnitude. We formalize 
this limitation as follows. Let 


5 


1 Max 

2 (f , 77 ) 6 2 




(A.21) 
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Then we require 


6 << 1 


(A.22) 


Define this diameter d of S by 


d = 2 


Max 
V) e 2 




and the height cj of S by 


Max 

(I T7)eZ 




It can be shown that 8 is an upper bound for co/d; hence, equation (A.22) implies that the 
ratio of the height of S to its diameter is small. Equation (A. 20) can always be ensured by 
a sufficiently fine grid. As a practical matter, we have adopted 

6< .066 (A.23) 


as a “rule of thumb” governing panel density. For a two dimensional cylinder, equation 
(A.23) would imply a maximum of 30° subtended angle per panel or a minimum of 12 total 
panels. 

Although it is not sufficient, equation (A.22) helps to guarantee that S is a close approxima- 
tion to the true surface. It also allows the expansion and subsequent integration in closed 
form of the induced velocity kernels. Moreover, it is doubtful that our assumed linear or 
quadratic distributions of singularity strength on S would be adequate for more highly 
curved panels. 
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APPENDIX B 

SINGULARITY STRENGTH DEFINITION 

B.l INTRODUCTION 

In section 4.2 it was noted that equation (13) represents an explicit solution to any of the 
boundary value problems described in section 4.1, once the source strength distribution 
a(0) and doublet strength distribution m(Q) are determined such that the specified boundary 
conditions are satisfied. 

In this appendix we shall describe how an approximation to the true network singularity 
distribution is defined for the purpose of achieving a numerical solution to the flow problem. 
The true singularity distribution will be approximated by a truncated Taylor’s series on each 
panel. Such a representation is valid on any interior part of the network providing the 
paneling there is sufficiently fine. It is less valid at a network edge where the true singularity 
distribution may become unbounded, in which case some error on panels adjacent to the edge 
will be unavoidable regardless of panel density. The approximation here as well as in the 
network interior can be improved by using a higher order distribution than currently 
employed by first order methods (which use constant source and doublet distributions). 

The next logical step as far as sources are concerned is to use a linearly varying source 
distribution on each panel. The order of doublet distribution consistent with such a source 
distribution is quadratic. This can be seen from equation (C.IO) which shows that it is the 
gradient of doublet strength which performs a function similar to that of a source, i.e., it 
creates a jump in a certain tangential component of velocity across a singularity surface. 

Hence, the gradient of doublet strength should also be linear, following that the doublet 
strength itself must be quadratic. 

In this report, we will consider only a linear singularity distribution on source panels and a 
quadratic distribution on doublet panels. There may be an advantage in using even higher 
order distributions, but as pointed out in reference 13, it would then be necessary to 
consider a higher order panel geometry definition for the sake of consistency. Specifically, 
we assume that the singularity strength X at a point (f , 17, f) on a panel S is given by 

Source : X(^, ??, f) = a(^ t?) = J t? (B. I ) 

Doublet: X(^, t?, f t?) 

= Mo + M j 7 ? + J ir}+ \ r?2 (B.2) 


Here (J, 7 ?, f) are local panel coordinates as defined by equation (A. 19). 

The unknown source and doublet coefficients (which may also be referred to as degrees of 
freedom) on the right sides of equations (B.l) and (B.2) respectively, are not assumed 
independent; rather they are linear combinations of an independent set of singularity 
parameters X| , X 2 , X^, . . . values are to be determined from application of 

the boundary conditions. The complete set of these singularity parameters will be denoted 
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as A. In otherwords, for each panel S we have 

C = BX _ (B.3) 

where C is the column vector of either the source or doublet coefficients and X is a column 
vector whose elements Xj^ form a subset of A. . In the next s^tion we describe; 1) how the 
set A is chosen for each type of network, 2) how the vector X is chosen for a panel, and 
3) how the rectangular matrix B is determined. (For a source panel B has dimensions 
3 X Ng; for a doublet panel B has dimensions 6 x Np. Ng and Nq are the size of X for 
source and doublet panels respectively and depend on the network type and on the location 
of the panel within the network.) 

B.2 CALCULATION OF DISTRIBUTION COEFFICIENTS 

The singularity parameters A are chosen to be the source or doublet singularity strengths 
X(f, Tj, f) at a set of discrete points on the network surface. The choice of S2 for the 
various types of networks employed is shown schematically in figure B.l. 

The circles represent points in J2. The intersections of the lines correspond to grid points 
of the network, i.e., panel corner points (see fig. A.l). A circle at one of these intersections 
therefore represents a grid point. A circle midway between two adjacent intersections 
represents the average of the corresponding two grid points. Finally, a circle centrally 
located among four adjacent intersections repjesejits the average of the four corresponding 
grid points. Here the average P* of N points Pj , P 2 , . . . Pisi is defined by 

= (P j + P 2 + . . +Pn)/N (B.4) 

and this definition holds even if some of the are identical. Technically, some of the points 
of S2 may not lie on any panel of the network since in general the panels defined in 
appendix A may not contain the grid comer points and/or midpoints; however, this is of 
no consequence since only the projections of the points in onto planes tangent to the 
panel surfaces are employed in future computations. The reasons for the choice of singularity 
parameter locations for various networks will be discussed in connection with control point 
locations in appendix C.4. We also note here that even numbered networks are doublet 
networks and odd numbered networks are source networks. There are no type 7 and 9 
networks. 

The singularity parameters in are ordered using the row index M = 1, 2, ... I as an inner 
index and column index N = 1 , 2, ... J as an outer index in a manner similar to a doubly 
indexed dimensioned variable in FORTRAN. The singularity parameters in A corresponding 
to two or more points of SI which physically coincide due to triangular panels, are identified 
in the ordering. This device can be used even when points do not coincide to create 
networks with fewer degrees of freedom than those in figure B.l. For example, the two 
wake networks are obtained in this manner from the type 6 (doublet/ analysis) network. 

The first is designated type 8 (doublet/wake #1) and is obtained by identifying all parameters 
in each column with the parameter at the head of the column. The second is designated 
type 10 (doublet/wake No. 2) and is obtained by identifying all parameters with the first 
parameter. 
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Type 1 (Source/analysis) Type 2 (Doublet/analysis) 



Type 5 {Source/design No. 2) Type 6 (Doublet/design No. 2) 



Type 8 (Doublet/wake No. 1 ) Type 10 (Doublet/wake No. 2) 

Figure B. 1. - Schematic of Singularity Parameter Location 
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Wejiow define the vector X of equation (B.3) for an arbitrary panel S. The components 
of X are the parameters Xj^e A corresponding to points Pj^ which belong to S or else 
are adjacent to points belonging to S. For this purpose a point P e is considered to belong 
to a panel if it is the average point of one or more of the panel’s comer points. (Hence some 
points in belong to more than one panel.) Also, a point Q e is considered adjacent to 
a point P belonging to S if each of its I and J indices differ by one at the most. 

The matrix B is obtained as the result of fitting the distribution of equation (B.l) or 
equation (B.2) to the singularity parameters of X by the method of least squares. To be 
specific we minimize the quantity 

5 (B-5) 

k 

with respect to the degrees of freedom in equation (B.l) or (B.2). In equation (B.5) the 
sum is over the N components of Xj^ which belong to S, and (f f j^) are the 
coordinates of the corresponding points Pj^ e expressed in the local coordinate system of 
the panel S. The weights are chosen to be unity unless ^ happens to belong to S 
in which case Wj^ is chosen large (nominally 10^). If we let C be th^vector having the 
coefficients (degrees of freedom) of X(J , 17, f) as its components and Vj^ the vector with 
components (1 , t?j^) in the case of a source distribution or 

2 ^k^» ^k^k» 2 

in the case of a doublet distribution, then the quantity R is minimized with respect to 
the coefficients when 


C = 


2 

k 


»k Vk V 


2 WkXkVk 
k 


Comparison with equation (B.3) reveals that 


B = Bf’ B2 


(B.6) 


(B.7) 


where 


B,. L 
k 


W V V T’ 
k k k 


(B.8) 


and Bn is the N column matrix whose kth column is 

Thus, for each panel S, the coefficients in equation (B.l) and (B.2) are expressed in terms 
of the Xj^ belonging to the panel. 
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B.3 CONTINUOUS DOUBLET/ ANALYSIS DISTRIBUTION 


The network singularity distributions defined in the previous section are basically splines 
although they lack the continuity characteristics usually associated with splines. Ideally, 
source strength and doublet strength and gradient should be continuous across panel edges. 
This is because discontinuities in source strength and doublet gradient across a panel edge 
induce logarithmically unbounded velocities there, and a discontinuity in doublet strength 
induces a jump in potential at the edge as well as a velocity which becomes unbounded as 
the inverse of the distance to the edge (see app. D.5). The problem with such flow 
anomalies is that flow boundary conditions imposed at control points close to these 
anomalies will concentrate on eliminating the anomalies in order to produce finite flow. 

The only way the flow anomalies can be eliminated is by greater continuity of the source 
and doublet splines ; hence, the boundary conditions will interact with the splines rather 
than control the finite flow in an appropriate manner. The weak flow anomalies produced 
by source and doublet gradient discontinuities are of little consequence to boundary 
conditions applied at panel centers; however, the strong flow anomalies produced by 
discontinuities in doublet strength can be of concern where panel center control points lie 
close to panel edges, as might be the case when the upper surface of a thin wing has 
different panel spacing than the lower surface. 

The least square procedure employed in the previous section does produce virtual continuity 
of doublet strength in regions where paneling is sufficient. In regions such as thin wing 
leading edges where the singularity strength gradients are large and the paneling is usually 
too coarse for the quadratic approximations to hold discontinuities do arise. It must be 
remembered that discontinuities in doublet strength correspond to jumps in potential 
which are not reflected in the gradient of the potential. Hence, calculation of velocities 
from the gradient of the doublet strength, as in example 4 of section 5.4.2, as well as the 
calculation of ACp for infinitely thin wings in many examples of section 5, can be erroneous 
without taking into account the fact that the variation of potential may be “lumped” at 
panel edges. This can be done after solution, i.e., after the singularity parameters are known. 
It is then possible to average doublet values along common panel edges creating a unique 
definition of doublet strength along grid lines. The distribution of doublet strength on any 
panel may then be modified by fitting the distribution coefficients to the values on the panel 
perimeter in a weighted least square sense. This refinement has been implemented in the 
pilot code and is responsible for more accurate values of Cp in cases involving the least 
square doublet/analysis spline of the previous section where panel density is sparse, 
however, in view of the preceding paragraph it would be far more desirable to construct a 
doublet spline with inherent continuity across panel edges. 

It is virtually impossible to construct a quadratically accurate doublet spline with exact 
continuity across all panel edges of a network when the distribution on a panel has only 
six degrees of freedom (equation (B.2)). The best we can do is achieve continuity at 
certain points along panel edges, e.g., at comer points. The type 4 (doublet/design no. 1) 
network is continuous at panel corner points and does a good job on the inverse problem, 
i.e., obtaining potential from gradient data. 
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In this section we will develop an alternate spline for the type 2 (doublet/analysis) network 
with sufficient continuity to yield more accurate velocity data at panel centers in cases 
where panel density is sparse. 

There are many ways to construct a doublet /analysis spline which is continuous at corner 
points. However, we wish to preserve certain desirable characteristics of the type 2 least 
square spline of the previous section. In particular we want the quadratic distribution on a 
panel to depend linearly on local singularity parameters only, on the nine parameters in an 
immediate neighborhood. Continuity at corner points then requires that a comer point 
value of the spline depend on the four adjacent singularity parameters only. This eliminates 
the choice of point values of the spline for singularity parameters as in the previous section, 
since four point values cannot in general be interpolated to obtain a fifth by a second order 
accurate formula. On the other hand, the following choice of singularity parameters does 
allow the second order accurate detemiination of corner point values by four adjacent 
parameters in the case where grid lines are straight. On any panel S we define the associated 
panel center singularity parameter in terms of the local doublet distribution by the formula 

4 

x=^ S (Pi + 2-Pi)] 

i= 1 

Here 

fi= Pift.ii.fi) 

is the ith panel corner point as shown in figure (B.2) with P 5 = P| and = P 2 . The 
distribution /i(f, is defined by equation (B.2) and 77 ) is the surface gradient of /i 
defined by the formula 

Tj) = T?) , Tj) , 0 J (B. 10) 

where 

M^(5. 7?) = + /^ jtj’J 

and 

Any network edge singularity parameter can be defined by equation (B.9) as well by 
allowing S to collapse to the panel edge or corner associated with the parameter 
(see fig. B. 1 ). For example, if ?] is a network comer point, then the appropriate 
expression for the singularity parameter X at?^j is obtained by sending ? 2 > ^3 and ?4 to 
P| ; hence, from equation (B.9) 

XatPj =p({j,t7i) (B.ll) 
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Figure B, 2 . — Pane! Schematic for Singularity Parameter Definition 


If P| and ?2 are network edge grid points, then the appropriate expression for the ^ 
singul^ity parameter X at the midpoint of and ?2 obtained by sending ^3 to ?2 
P4 to Pj ; hence, from equation (B.9) 

Xat i (?, +^ 2 ) = T )+ • p2-^l)] 

+ I [m(^, 7?2) + ^ v/i(f2,r?2) . (P 1 -P 2 )] (B.12) 

Since 17) is quadratic along the line joining Pj and?2 expressions in square brackets 
on the right side of equation (B.l 2) can be shown to be equal so that 

Xat j(p,+P2)=M(S],r„) + • f 2 -Pl) 

= (b’ ^2) + ^ (^2> ^2) * (^1 - P2) (B.l 3) 

The last expression is closely associated with the definition of the singularity parameters 
used to solve the one-dimensional quadratic interpolation problem. Assume we wish to 
interpolate values M , ^n + 1 points k = 0, 1 , . . . n + 1 where 

^0 = ^1=2 (""o + ’^l) > ^2 = i (xi + X 2 ) = x„ 

as shown in figure B.3, It is possible to define a function /i(x) such that 1) ii(x) is quadratic 
on each interval [xj_ j, Xjj, i = 1, 2, . . . n; 2) jL/(x) is continuous and 

g '(x) = is continuous on the full domain [x^, x^^J ; and 3) g(x) interpolates the 

values = 0, 1, . . , n + 1. Moreover, the g(x) with such properties is unique. 
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Figure B.3. — One Dimensional Quadratic Interpolation 
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Equation (B.I 6 ) has the following interpretation; plot the in place of the and connect 
the values with straight lines. Then and M at are the value and slope at xj^ of the 
straight line segment connecting Xj, and Xj^.j. j. From equation (B.16) we can obtain the 
coefficients of the quadratic doublet distribution in each interval in terms of the Xt. We 
have 

Mk(x) = Mk + Mk (x-?i)+^Mk (B.17) 

where 



and 

dj = (fk-fk-l)- d2 = (xk-Xk-i). d3 = (fk+l-^k)- 

Note that the coefficients depend only on the singularity parameters belonging to the interval 
and its immediately adjacent intervals. The Xk can be obtained in terms of the Mk by solving 
an equation set consisting of Xq = Mq. + 1 ^ ^n + 1 *be first equation of (B. 1 8) for 
every interval, using a tridiagonal equation solution algorithm. 

The definition of equation (B.9) is an attempt to generalize equation (B.15) for the case 
of a surface spline. As in the one dimensional case, the next task is to determine the comer 
point values of doublet strength and gradient. Given any comer point ?(I, J) of the network, 
we define a local (^, tj, f) coordinate system with origin at the corner point and with the 
(J, 77 ) plane approximately tangent to the true network surface. The f, T unit vectors are 
defined as follows: 



A A (P3I -^ll) (^32-^12) (^33-^13) 

*|(p31-Pii)+2(?32-?,2)+(?33-^3)| 


(B.19) 



A A A 
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Here P'»2 •« the point P(I, J) and the other Py are the eight adjacent corner points defined by 

P,,=P(I-1,J-1) 

P2, =P(I, J - 1) 

?3j = P(I + 1, J - 1) 

P,2 = P(I- 1, J) 

P32 = P(I + 1, J) 

Pj3 = P(l- 1, J + 1) 

P23 = P(I, J+ 1) 

P33 = P(I + I, J + 1) (B.20) 

(If any index exceeds a grid point row or column limit it is assumed to be replaced by 
that limit.) These comer points are shown schematically in figure (B.4) along with 
the four relevant singularity parameters. 



Figure B.4. - Corner Point Schematic for Spline Corner Point Value 

Determination at the Corner Point P22 Common to Four Panels 
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Let and be the doublet strength and gradient at the origin ? 22 - These quantities are 
then to be determined in terms of the four Xj, by minimizing the function 


• (Sj + Vl +dj) 

1 * 


with respect to 




Here 


^ 1 " ^32 ' ^ 22 ’ 



(B.22) 


and 





= d 


1 


and 

(ii' 1i.fi)- Si 

The weight W is nominally chosen as 10^. The weight is chosen as the square of the 
5 component of the^ and with largest magnitude, the weight is chosei^imilarly. 
The choice of R in equation (B.2 1) is motivated by an attempt to obtain 
fitting a quadratic distribution to Xj, X 2 , X^, X^ with aj, a 2 , aj being the second order 
coefficients. Such a fit has two degrees of freedom, however, if is parallel to^^ and^2 
is parallel to^ 4 , these degrees of freedom affect the second order coefficients only and 
the values of and^p^. are determined uniquely by the Xj. 

The function R is minimized by values of VMc> ^2’ ^3 "'hich can be obtained from 
a formula similar to equation (B.6), i.e.. 
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"V 


r^c 1 



32 


^3 


Known Matrix 
Determined Solely From 
Geometry of Py 



(B.23) 


As in the one-dimensional case, once the dependence of all comer point values and gradients 
on the network singularity parameters has been found the panel distributions 77 ) may 
be calculated. Return to figure (B.2) and equation (B.2). The six coefficients of v) in 
equation (B.2) are determined as follows. We require that the distribution actually attain 
the corner point values just computed, i.e,, 

M(?i, T?i)=Mi ; i=l,...4 (B.24) 

where ;xj is the value of at we also require that equation (B.9) hold for the panel 
center singularity parameter and any network edge singularity parameter associated 
with the panel. All remaining degrees of freedom (if any) are then eliminated by 
requiring that 

VM(ii,T?i)= VMi ; i=l,...4 (B.25) 

hold in a least square sense, where is the value at^j. Specifically, for a panel 

in the network interior, the six coefficients of t?) are obtained by minimizing the 
function 


+ ■ (Pi«-P|)]^ (B.26) 

with respect to 

Mo> 
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where and are the ^ and t? components of V/ij, and X is the panel center singularity 
parameter. The weight W is nominally chosen as 10^. The weights and are chosen as 
the maximum values of |^jl and lT?j| respectively. Only the panel center singularity parameter 
defining equation is displayed on the right side, but for network edge panels the associated 
network edge singularity parameter defining equations should also be included. Note that 
for network corner panels the doublet distribution can only achieve specified corner 
point values and singularity parameters in a least square sense. The procedure for minimizing 
R is similar to that used to derive equation (B.6) and the result is a B matrix to replace 
that of equation (B.3). The vector C is again the coefficients of equation (B.2) and the 
vector X consists of the singularity parameters which belong to S or immediately adjacent 
panels. The 6 x Nj) matrix B again depends only on the geometry of the corner points of 
S and those of its immediately adjacent panels. 
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APPENDIX C 

CONTROL POINT DEFINITION AND 
BOUNDARY CONDITION SPECIFICATION 


In this appendix we describe the selection of a set F of control points on a network. (By 
definition control points are points at which boundary conditions are applied, i.e., points at 
which the integral equations or auxiliary conditions of the problem hold exactly.) In 
addition, we describe the form in which the boundary conditions may be specified. 

C.l CONTROL POINT LOCATION 

The set F for each network is shown schematically in figure C.l . (There are no type 7 and 
9 networks.) 

The circles represent points in F , As in figure B, 1 the intersections of the lines correspond 
to grid points of the network and the squares correspond to panels. Circles located on 
edges or at intersections correspond to the same network locations as in figure B. 1 —at least 
for the present. However, a circle in the middle of a square (which we call a panel center 
control point) denotes the origin of the local coordinate system (see app. A) of the 
corresponding panel rather than the average of the four adjacent comer points. The points 
in F may be ordered using the row index as an inner index and the column index as an 
outer index. For this purpose points which coincide with previous points are deleted. 

Hence, F and A have the same size for network types where A and are schematically 
similar, i.e., types 1, 2, 8, and 10. This may not be the case with the design networks when 
triangular panels are present. We do not specifically exclude triangular panels from design 
networks, but we do require that the presence of triangular panels result in sets F and A 
of identical size. For example, the collapse of an entire edge of a type 4 network would be 
permissible. 

Once F has been ordered, it is necessary to withdraw the control points on the network 
edges slightly to avoid numerical difficulties associated with infinite self-induced velocities 
(see app. D.5). The displacement is accomplished in the following manner. First, assume 
the control point P is one of the four corner grid points of the network. Without loss^ 
of generality, we can assume P = P( 1 , 1 ) (see fig. A. 1 ). Let N be the unit normal at P of 
the panel containing P. The control point P is then redefined by the formula 

P = P(l, D + a^Dj + 02 ! (C.l) 


where 


Dj = N®[P(1,2)-P(1, 1)] 
D2 = N® [P(l, l)-?(2, 1)] 

d,-d,/15|| 

Dj - Dj/ISjI 

« = e(|D,|.|5j|) 
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Figure C. /. — Schematic of Control Point Locations 
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Here e is a small number (nominally 1/4 x 10“ ^). If either P(l, 2) or P(2, 1) coincides 
with P( 1 , 1 ) we replace it by P(2, 2). 

Next assume the control point "P is located at the edge of the network midway between two 
^rid points. Without loss of generality we can^assume ?= y (P(l, 1) + ?(1, 2)). Again, let 
N be the unit normal of the panel containing P. The contra point P is then redefined by 
the formula 


-► 1 r-> -*■ 1 ^ 

P=y [P(l, l) + P(l,2)J+aDj 


(C.2) 


where 

Dj = N ® [P(l,l)-P(l,2)] 

D|. D, /|d,| 
a = € |d J 

Here e is a small number (nominally 10“^). If P(l, 1) coincides with P(l, 2) the displace- 
ment procedure of equation (C. 1) is used. 

Note that the displaced network edge control points do not necessarily lie on their respective 
panel surfaces. However, this causes no problems since induced velocities at these points 
are computed in a special way taking into account only edge induced effects. (See the 
discussion in equations (C. 1 6) and (C. 1 8).) 

C.2 PANEL CENTER BOUNDARY CONDITIONS 


In this section we describe the form of the boundary condition equations at panel center 
control points. We first consider boundary conditions involving velocities where the 
equations may take the form 


* ^ l ) ~^n 


(specify cy, cy, Pn) 


(C.3) 


(% ■ ^u) + (\ • Vl) = (specifyTy, 


(C.4) 


Here Vy is the total upper surface velocity (perturbation plus free stream) and Vy is the 
total lower surface velocity. By upper surface we mean the side of the panel surface on 
which lies the positive f axis of the local coordinate system. The lower surface is the 
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other side. The unit vector n is the upper surface normal and ty and ty are surface tangent 
vectors. (Arbitrary vectors may be specified but they are projected onto the tangent plane of 
the panel at the control point before application.) Conditions of type (C.3) are employed on 
analysis networks while conditions of type (C.4) are used on design networks. These condi- 
tions cover problems involving both interior and exterior flows as well as thin sheets. 

To explain the application of the equations (C.4) and (C.3) we first consider a panel of the 
boundary surface B which bounds D on one side only (e.g., the wing section shown in 
fig. 21 ). Without loss of generality we can assume this side is the upper side. Lower surface 
^^locities are^hen irrelevant and we may assume Cy = 1 , Cy = 0 in equation (C.3) and 
|ty| = 1 and ty = 0 in equation (C.4). Analysis boundary conditions then reduce to the 
specification of the fluid velocity normal to the panel at the control point. For example, 
conditions of type (C.3) with cy = 1 , cy = 0, = 0 are the usual analysis boundary condi- 

tions on impermeable surfaces. On the other hand, both the unit tangent vector T and the 
^elocity component in the direction ty must be specified for design conditions. Often 
ty is selected to be the stream direction at the control point obtained from analyzing the 
flow about the existing geometry and |3^ is the desired velocity magnitude to be produced 
by a perturbed geometry. 

Next, let us consider a portion ofithe boundary surface which bounds D on both sides 
(e.g., the model of fig. 19). In order to control the flow on both sides, the surface must 
conceptually be represented by superimposed source and doublet networks. Without loss 
of generality we can assume the orientation of both networks is the same regarding upper 
and lower surface designations. We also assume that the paneling for both networks is 
identical so that a pair of boundary conditions is applied at each center control point 
location. This pair may consist of two analysis, one analysis, and one design or two design 
type boundary conditions. To study these combinations it is convenient to define a set of 
boundary conditions equivalent to (C.4) and (C.3). For this purpose we define an “average” 
velocity by 



and a “difference” velocity by 


V 


D 



(C.6) 


Then (C.4) and (C.3) are equivalent to 


^specify c^, Cj), j (C.7) 


Fa ' ) + (to ■ Vd ) = (specify t^, tp, ) 


(C.8) 
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respectively, where 




Equivalence is established by noting the inverse transformation, i.e., 


Vu = 

-► ”► 1 “► 
Vl = 


1 

1 

^L“ 2 

A 

-¥ 1 ^ ^ 
*L ■ 2 


(C.9) 


Let us first consider the case where both boundary conditions are of analysis type, i.e., of 
the form (C.7). If the boundary value problem is well posed, the boundary conditions must 
be independent and consistent, hence, without loss of generality we can assume that one 
condition controls the normal component of the average velocity (i.e., cp = 0) and the 
other controls the normal component of the difference velocity (i.e., c^ = 0). In theory, 
it doesn’t matter which boundary condition is associated with the source panel control 
point or doublet panel control point. However, in the present numerical development 
(i.e., in the existing pilot code logic) a distinction must be made. The present developmejit 
calculates the influence coefficients for directly (see app. D.3) but obtains those for Vp 
by using the formula 


A 

Vp = an + v#x 


(CIO) 


which follows from equations (15), (16), and (B.IO). 

Here o is source strength and the doublet strength gradient at the control point. 

Because of the logical structure of the present method (in which network independence is 
maintained) only the first term on the right side of equation (C. 10) is considered for a source 
panel control point and only the second term is considered for a doublet panel control 
point. It is then essential that the boundary condition controlling the normal component 
of the difference velocity (i.e., c^ = 0) be associated with the source panel control point 
since only for a source panel control point is this component computed correctly. (Note 
n = 0, see equation (B.IO).) 
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A similar situation exists in the case of one analysis boundary condition (type (C.7)) and one 
design boundary condition (type (C.8)). Here the analysis boundary condition must be 
associated with the source panel control point and the design boundary condition with the 
doublet panel control point. The case of two design boundary conditions is slightly different 
since it is no longer possible to assume withouHoss of generality that one boundary 
condition controls a tangential component of and the other a tangential component of 
Vj). Such an assumption is not too restrictive since one can probably ^chieve ^y general 
design result of this type by controlling the tangential components of and Vq separately. 
At any rate, it is clear that^nder such an assumption the source panel boundary condition 
should be used to control and the doublet panel boundary condition should be used to 
control Vq. 

As a final note, equation (C. 10) implies that a surface across which normal velocity is con- 
tinuous may be represented by a doublet network alone. On the other hand, a surface 
across which tangential velocity is continuous may be represented by a source network 
together with a doublet network of constant strength. The latter network is unnecessary 
if there is no jump in potential across the surface. 

Next, we consider specification of potential at panel center control points where the boundary 
condition equation may take the form 

ku</>u + ^specify ky, kL,^^) (C.ll) 

Here and upper and lower surface values of perturbation potential respectively. 

Boundary conditions imposed by equation (C. 1 1 ) are design (Dirichlet) conditions, but 
differ from those imposed by equation (C.4) in that scalars rather than components of vectors 
are specified. For this reason the symmetrically defined singularity parameter and control 
point locations of network types 1 and 2 are most suitable for applying these conditions. To 
avoid constructing additional design networks with identical properties we simply use these 
network types substituting equation (C. 1 1 ) for (C.4). 

The application of equation (C. 1 1 ) is obvious from our earlier discussion. We note only the 
analogues to equations (C.5), (C.6), (C.7), (C.9), and (C.IO). We have 

A~ ^ O’ (C.12) 


4>\j-4>A 




7^D' 


“ 2 


(C.13) 


and 


kA^A ^ (specify k^, kp, 


(C.14) 
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where 


1 

1 

‘^L=2^A’*'D 


Finally, we have from equation (15) 

= ^ (C.15) 


C.3 NETWORK EDGE BOUNDARY CONDITIONS 

Boundary condition equations of type (C.7) with = 1 , Cp = 0 and 0 are automatically 
imposed at all edge control points of the doublet/analysis network. Because of the singular 
behavior of the velocity induced by a finite doublet distribution at a panel edge, this boundary 
condition in effect controls the continuity properties of the distribution across the edge. To 
give an example we note (from equation (D.141)) that near the common edge of two smoothly 
adjoining doublet networks the downwash (normal velocity) is given by 

(Va • n) = ^ + 2 Am log(|7?l) + regular terms j (C.16) 

where t/ is the tangential coordinate perpendicular to the edge, Am is the Jump in doublet 
strength across the edge and Am^ is the jump in the derivative of doublet strength in the 
direction perpendicular to the edge. A control point placed near the edge requiring that 
downwash be finite will tend to make Am vanish (because of the strong antisymmetric 

singularity — ), i.e., m continuous across the edge. A similar control point on the opposing 
V 

panel of the adjoining network will, in addition, force Am to vanish (because of the weaker 
symmetric singularity logdr?!)), thereby establishing continuity of M- 

Because of small, unintended discontinuities in geometry between panels of adjacent net- 
works due to the approximate nature o^thc surface fit technique, such singularity matching 
cannot be accomplished dependably if is computed in an exact manner. (See section 
6.3.5 for the effect of geometry gaps on doublet strength matching.) Consequently, a 
special algorithm is used. This algorithm considers only velocities induced by panel edges 
adjacent to the control point. Moreover, the algorithm computes an edge induced velocity 
as if the edge were the straight line joining the corresponding two grid points rather than the 
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actual panel edge. This means that two doublet network edges are considered physically 
joined from a potential flow standpoint if and only if the rectilinear curves formed by connecting 
their edge grid points coincide. Unintended discontinuities in surface slope between adjacent 
panels of two abutting networks can also occur because of the approximate nature of the 
surface fit algorithm. However, the effects of such an error remain quite local in character and 
no special modification is required. 

The edge boundary conditions for the remaining doublet networks are the same as for the 
doublet/analysis network and they control continuity of ^ and/or m ' as well. As discussed 
in section 4.1 these edge conditions are fundamental to the boundary value problem for design 
networks and fulfill the auxiliary conditions of type equation (12). Network type 4 assumes 
inflow on the left and lower edges whereas network type 6 assumes inflow on the lower edge 
and outflow on the upper edge. The wake network types 8 and 10 are special design networks 
used in place of a regular doublet/design network when a reasonable guess of the direction of 
velocity is deemed sufficient. (A common use is for the representation of vortex wakes.) 

Auxiliary conditions corresponding to equation (11) can be substituted for the continuity 
boundary conditions at the edge control point heading each panel column in networks type 4 
and 6. These conditions take the form 

^cu(n • Vu) + CL (n • Vl ) dS = 0 (C. 1 7) 

C 

where the integration extends over a panel column and Cy and Cl are specified. The integral 
is evaluated by summing over each panel in the column the product of the panel area and 
the integrand evaluated at the panel center control point. The networks are assumed to be 
paneled with columns aligned along streamlines. In practice this requirement may be relaxed 
considerably, and only a general streamwise alignment appears necessary. Technically, new 
network types should be assigned for application of these closure conditions; but to avoid 
constructing additional design network with identical properties we simply use types 4 and 
6 with the above mentioned modification. 

The edge boundary conditions for the source/design networks 3 and 5 perform the same 
functions as those for 4 and 6. For application of equation ( 1 2) a different numerical device 
is used which takes advantage of the fact that a finite strength source panel induces an 
infinite tangential component of velocity at its edges. To give an example we note from 
equation (D. 121) that near the edge of two smoothly adjoining source networks the tangential 
component of velocity perpendicular to the edge is given by 

(Va ■ t ) = ^ 2 ao log (|r?|) + regular terms) (C.I8) 

where t is the unit tangent vector perpendicular to the edge, rj is the coordinate in this 
direction measured from the edge and Aa is the jump in source strength across the edge. A 
control point placed near the edge requiring that (Va • ^ be finite will tend to make Ao 
vanish, i.e., source strength continuous across the edge. This in turn will accomplish 
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equation (12). The same special algorithms for evaluating for doublet edge control 
points are used here. Finally, equation (C.17) is applied in the identical manner for network 
types 3 and 5 as for network types 4 and 6 when closure is preferred. 

In the case where specification of potential is substituted for specification of normal velocity 
at center control points of network type 2, potential is also specified at all edge control 
points. 

C.4 RATIONALE FOR CONTROL POINT/SINGULARITY PARAMETER LOCATIONS 

The set ^2 of singularity parameter locations and the set F of control point locations for 
various network types are shown schematically in figures B.l and C.l respectively. Let us 
first consider the type 1 (source/analysis) network. We note the following: 1) The sets 
r and f2 are schematically the same for this network (although the precise point locations 
are slightly different); 2) The points in F (or S2) are symmetrically located with respect 
to the network grid point schematic; and 3) The points in F and H are located at the 
panel centers. 

The locating of the points of F at panel centers is quite reasonable in view of the discussion 
at the beginning of section B.3. Both the panel geometry and singularity splines are analytic 
at panel centers, and moreover, the panel centers are the points which are farthest from the 
flow anomalies induced by discontinuities in geometry and singularity strength at panel edges. 
The locating of the points of at panel centers is acceptable from two points of view. 

First, ^2 and F have the same size leading to the same number of equations as unknowns. 
Secondly, there are a sufficient number of local singularity parameters (i.e., singularity 
parameters belonging to a given panel and its immediately adjacent neighbors) to determine 
the three coefficients of the panel source distribution via the method of least squares. The 
smallest set of local singularity parameters occurs for each of the comer panels in a net- 
work and the four local parameters in the set are sufficient to determine three coefficients. 

The fact that the points of F and J2 are symmetrically located is reasonable considering 
the allowable boundary condition types for this network. These boundary conditions can 

involve ^ only, and hence, involve no preferred direction on the network surface as 
opposed to design boundary conditions which involve directional derivatives of 0. 

The fact that the points of F and J2 are identical is among other things due to considera- 
tions of spline stability. As with any other application of splines, we desire that our source 
and doublet splines be stable relative to the type of boundary conditions applied. This 
means that small changes in the boundary conditions should result in similarly small changes 
in the singularity parameters. For given sets F and J2 the spline may be stable with respect 
to one type of boundary condition but not another. Also, for a given set ^2 and a given type 
of boundary condition, the spline may be stable for certain control point locations but not 
others. A discussion of spline stability is beyond the scope of this report (see e.g., ref. 27); 
however, the source/analysis spline can be shown to be stable with respect to application 

of boundary conditions at panel centers, using equation (C. 10). From experience, it also 
dn 

appears to be stable with respect to the application of 0 boundary conditions at panel centers. 
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The rationale for locating the points of the sets V and in the case of the (doublet/ 
analysis) type 2 network is the same as for the (source/analysis) network. The addition of 
network edge singularity parameters is due primarily to the quadratic order of the doublet 
distribution on each panel. We note that the six coefficients of such a distribution 
could not be determined from the four local singularity parameters in the case of a network 
corner panel, were the source/analysis spline to be used. With the singularity parameter 
arrangement shown in figure (B,l) for the doublet/analysis network type, the local set of 
singularity parameters for every panel contains nine singularity parameters which is more 
than adequate to determine six distribution coefficients by the method of weighted least 
squares. The corresponding addition of edge control points to P keeps the number of control 
points equal to the number of singularity parameters and allows for doublet strength and 
gradient matching across network edges. It can be shown using equation (C.15) that the 
arrangement of points in P and 12 in the case of the doublet/analysis network type is stable 
relative to the application of 0 boundary conditions. From experience, it also appears 

A/A 

stable relative to the application of — type boundary conditions on a surface bounded by 
the fluid domain on both sides. 

The control points for network types 3, 4, 5, and 6 are not symmetrically located due to 
the directional nature of design type boundary conditions. The discussion corresponding to 
figure 4 indicates the necessity of defining “auxiliary” boundary conditions at edge points 
of a design network where inflow is anticipated. For network types 3 and 4 we assume inflow 
across the bottom and left side edge (relative to the schematic displayed in fig. C. 1 ); for 
network type 5 we assume inflow across the bottom edge only; and for network type 6 we 
assume inflow across all but the top edges. The network interior control points are located 
at panel centers for the same reason as in the case of the analysis networks. Spline stability 
relative to design type boundary conditions then requires that the singularity parameter points 
be located at grid corner points in the case of the network types 3 and 4, and at edge midpoints 
in the case of the interior parameters for network types 5 and 6. In each case the number of 
singularity parameters equals the number of control points. 

The wake type networks 8 and 10 are in reality the design type network number 6 where the 
boundary condition has already been applied directly to the spline, resulting in fewer degrees 
of freedom. The type 8 network is obtained by applying to all type 6 control points, except 
for those along the bottom edge, the design type boundary condition that Vju be zero in the 
panel column direction. (Assuming that the panel column direction corresponds roughly to 
the freestream direction such a boundary condition implies that the ACp = 0 where Cp is 
computed using the linearized formula.) By integration we see that such a boundary condition 
implies that all singularity parameters in each column of the type 6 singularity spline will be 
equal to the value of the singularity parameter at the head of the column (i.e., at the bottom 
edge of the schematic). Performing such singularity parameter identification results in the 
type 8 network. The remaining control points and singularity parameters are used for doublet 
strength and derivative matching (see sec. C.3). The type 10 network is obtained by requiring 
in addition that Vfi be zero in the row direction along the first (bottom) row of control 
points, which can be achieved by identifying all singularity parameters in that row with the 
first. This results in a network where doublet strength is constant everywhere. Such a 
network creates a jump in potential but no jump in velocity (see equations (C. 10) and (C.15)). 
In contrast, the type 8 network creates a jump in the component of velocity perpendicular 
to the column direction. 
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APPENDIX D 

PANEL INFLUENCE COEFFICIENT GENERATION 

D.l INTRODUCTION 

in this appendix we shall calculate the potential and velocity induced by a source or doublet 
distribution on a curved panel. As shown in figure (D.l), let S be th^ curved panel surface, 

2 its tangent plane projection, Q a point on S,^ the normal to S at Q, and P a field point. 



The perturbation potential 0 at P induced by a source distribution of strength a on S 
is defined by 

s 


where 


R = (^-x,i?-y, f-z) = Q-P 


(D.2) 
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and 


R= |R| 




x)“ + (7?-y)^ + (f-z)^ 


The perturbation potential </> at ? induced by a doublet distribution of strength non S 
with doublet axis in the n direction is obtained by taking the directional derivative of 
in the Indirection: 



S 


The perturbation velocity v induced at P by a source or doublet distribution on 
defined by 


(D.3) 


S is 


V = V 0 


(D.4) 


We assume that the surface S is defined by equation (A.20): 

f = a?2 + (j, rj) e Z 


(D.5) 


We also assume that S does not deviate significantly from 2, more precisely that 


5 << 1 


(D.6) 


where 


6 


1 Max 

2 (I T?) e 2 


(7a¥ + bV 


(Nominally, we assume 6 < .066 see equation (A. 23).) 


(D.7) 


The distribution of singularity strength on S is assumed to be linear in the case of a 
source panel and quadratic in the case of a doublet panel. To be specific we assume 


a = OQ + ajf + a^rj, ({, j?) e 2 


(D.8) 
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and 


1 1 ? 

M = + + . (tTj)eS 


(D.9) 


which are the same as equations (B.l) and (B.2) respectively. 

D.2 EVALUATION OF SOURCE AND DOUBLET INTEGRALS 
FOR AN ARBITRARY FIELD POPSIT 

Let us first consider the evaluation of source potential defined by equation (D.l). Evaluation 
of source velocity and doublet potential and velocity will be quite similar. The first step 
in the evaluation procedure is to transfer the integral over S to the equivalent integral over 
2. (See discussion of fig. (A. 2) for a precise definition of the plane quadrilateral 2.) 

We have 

2 

From equation (D.5) we obtain 


A _ 1 

V 1 + 4a2f2 + 4b2rj2 


(-2a£,-2bT?, 1) 


(D.ll) 


and 


sec (t n) = v/l + 4a-^^ + 4b“Tj2 


(D.12) 


Substitution of equations (D.2), (D.8) and (D.12) into (D.IO) yields an explicit integral 
for 0. However, the integral cannot be evaluated in closed form as it stands. By employing 
the hypothesis that 6^ is negligible compared to unity (hypothesis (D.6)) the integrand can 
be approximated by terms that are integrable in closed form. A uniform approximation 
to sec (f, Tj) can be obtained by noting that 

4a-S-+4bV<l65- 


hence 

sec(f,n)«l (D.14) 
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A uniform approximation to 1/R is somewhat more difficult to obtain since this factor is 
singular. Let (xq, Yq, 0) be the point on S closest to (x, y, 0) (see fig. D.2) and set 



Figure D.2. - Geometry for Curvature Approximation 



Let 


Max I ■ ^ol ) 
a, r?)e2( r ) 


(D.17) 
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Then 


-2h(f-Zo) + (f-zo)- 


r- + h- 


2erh + 
r- + h“ 


< e + e-" 


(D.18) 


T . 


Therefore, if e~ is everywhere negligible equation (D. 16) yields 


1 1 h(r - Zq) ^ : 

-=s-+ , p= Jx~ + h‘ 


R P p3 


(D.19) 


But 


e = 


Max I yp)| 

(t, T?) e Z I r 


Max 

(f, 1?) c Z 


( + '‘o) (^ ' '‘o) + + yo) - yo)| 


Max 

(t, T?) e Z 


Ja2(t + xo)2 + b2(T) + yo)2 


J(t-xo)^+(’?-yo)' 


Max 

(^, T?) e Z 


|^a2(t + xo)^ + b2(7, + yo)2 . zj 


< 86 


(D.20) 


A much better bound on e is available when (x, y) is several panel diameters away from Z 
and the assumption that 6^ is negligible becomes unnecessary. However, in this case a far 
field expansion will be used to obtain an efficient approximation to the right side of 
equation (D.IO). Note that the panel curvature effect is contained in the second term on 
the right side of (D.19), this term being zero for a flat panel. 

Substituting equations (D.8), (D.14) and (D.19) into (D.IO) and rearranging them, we obtain, 
after considerable algebraic manipulation: 
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<t> = o(x, y) I( 1 , 1 ) + a^j(x, y) 1(2, 1 ) + Oy(x, y) l( 1 , 2) 


(D.21) 


where 

ct(x, y) = OQ + ajX + a^y 
a^(x, y) = 

Oy(x, y) = a^ (D.22) 

Note that a^(x, y) and Oy(x, y) are constant, although formal dependence on x and y 
has been displayed to emphasize the shift in the origin of the Taylor’s series expansion to 
the point (x, y). 

Here 


1(M, N) =~ I H(M, N, 1 ) + a[hH(M + 2, N, 3) + 2xhH(M + 1 , N, 3)] 

+ b[hH(M, N + 2, 3) + 2yhH(M, N + 1, 3)] 
+ c[hH(M,N,3)]j 


where 


^ ^ u 2 

c = ax- + by - Zr 


(D.23) 


(D.24) 


and 


H(M. N, K) = 


fj(^. 


(D.25) 


The H integrals will be evaluated in the next section. The leading term in the righthand 
side of equation (D.23) is due to the leading term on the righthand side of equation (D.19) 
and thus corresponds to a flat panel. The remaining terms having coefficients of a, b, and 
c are due to the second term on the right side of equation (D.19) and constitute panel 
curvature effects. 

To find *v, equation (D.21) can be differentiated. For this purpose Zq may be treated as 
constant, although formally zq depends on x and y because (xq, yg, 0) depends on x and 
y, being the point on 2 closest to (x, y, 0). However, the derivatives of zq with respect to 
X and y either cancel each other or are negligible because of hypothesis (D.6). The 
derivatives of the H integrals then are simple combinations of the H integrals themselves, i.e.. 
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— H(M, N,K) = -(M- 1)H(M- 1,N, K) + KH(M + 1,N, K+2) 
dx 


— H(M, N, K) = - (N - 1 ) H(M, N - 1 , K) + KH(M, N + 1 , K + 2) 

9y 


- H(M,N,K) = -KhH(M, N, K + 2) (D.26) 

dz 


It turns out to be easier to calculate v by differentiating equation (D.IO) to obtain 

— sec (f , n) dfdr? and using a generalized form of equation (D.19), 
^ 47tR3; 

that is 


rK 


1 Kh(f-ZQ) 


— + 




(D.27) 


Equation (D.27) can be obtained by raising equation (D.19) to the K power and expanding 
the righthand side by the binomial theorem. In either case we obtain 

V = ct(x, y) J(1, 1) + y) 1) + Oy(x, y) J(l, 2) (D.28) 


where 


J(M, N)= [jx(M, N), Jy(M, N), J^IM, N)] 


(D.29) 


i. e., V = o(x, y) Jjj(l, 1) + 
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and 


Jx(M, N) = - ^ j H(M + 1 , N, 3) + a[3hH(M + 3, N, 5) + 6xhH(M + 2, N, 5)] 

+ b[3hH(M + l,N + 2, 5) + 6yhH(M+ 1,N + 1,5)] 

+ c[3hH(M+ 1, N, 5)]j 

Jy(M, N) = -^ I H(M, N +1, 3) + a[3hH(M + 2, N + l,5) + 6xhH(M + 1,N + 1,5)] 

+ b[3hH(M, N + 3, 5) + 6yhH(M, N + 2, 5)] 

+ c[3hH(M, N +1, 5)]1 

J,(M. N) = - — I - hH(M, N, 3) + afH(M + 2, N, 3) - 3h“H(M + 2, N, 5) 

^ 47t I L 

+ 2xH(M + 1 , N, 3) - 6xh“H(M + 1 , N, 5)] 

+ b[H(M, N + 2, 3) - 3h-H(M, N + 2, 5) 

+ 2yH(M, N + 1 , 3) - byh^HCM, N + 1 , 5)] 

+ c[h(M, N, 3) - 3h-H(M, N, 5)] j 

Doublet potential and velocity can be evaluated similarly using equations (D.3), (D.l 1) and 
(D.27). 

However, the assumption (D. 14) is unnecessary since sec(f , t?) does not appear in the 
product ndS. We obtain 

4 > = p(x, y) 1(1, 1) + Mj^(x, y) 1(2, 1) + Py(x, y) 1(1, 2) 

y) 1(3, 1 ) + Pxy(x. y) 1(2, 2) + ^l^yyCx, y) K 1 , 3) (D.30) 
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where 


Mx, y) = Mo + MfX + M^y + ^ Mjjx 2 + ^ M^^y^ 

Mx(x, y) = M^ + Mf^x + Mj^y 

My(x,y) =M^ + M^Tj^+M^^y 

Mxx(x,y) =Mjj 
Mxy(x,y) =M^^ 

Myy(x. y) = Mtjtj (D.31) 

and 

I(M, N) =— |hH(M, N, 3) + a[H(M + 2, N, 3) + 3h2H(M + 2, N, 5)+6xh2H(M + 1, N, 5)1 
4ir I L 

+ b[H(M, N + 2, 3) + 3h2H(M, N + 2, 5) + byh^HCM, N + 1 , 5)] 

+ c[- H(M, N, 3) + 3h2R(M, N, 5)] j (D.32) 

As in equation (D.23), the terms containing a, b, and c as coefficients are the effects of 
panel curvature. 

We also obtain 

V = y) J( 1 , 1 ) + y) 1 ) + My(x, y) J( 1 , 2) 

+ ^ Mxx^’^' 1) + Mxy(x, y) J(2, 2) + ^ y^ J(», 3) (D.33) 

where 

t(M, N) = [j^CM, N), Jy(M, N), J^(M, N)] (D.34) 
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and 


Jx(M. N) 


Jy(M, N) 


Jz(M, N) 


= I 3hH(M + 1 , N, 5) + a[3H(M + 3. N. 5) - 2H(M + 1 , N, 3) 


+ 15h-H(M + 3, N. 7)- 2xH(M, N. 3) + 30 xh~H(M + 2, N, 7)J 
+ b|^3H(M + 1 . N + 2, 5) + 1 5h-H(M + 1 , N +2, 7) + 30 yh~H(M + 1 , N + 1 , 7)j 
+ c[-3H(M + 1,N, 5)+ 15h“H(M+ l.N. 7)]j 


= ^ |3hH(M. N + 1, 5) + a[3H(M + 2, N+ 1, 5) + 15h“H(M + 2, N+ 1, 7) 
+ 30xh-H(M+ 1,N+ 1,7)] + b[3H(M, N + 3, 5) - 2H(M, N + 1 , 3) 

+ 15h“H(M, N + 3, 7) - 2yH(M, N, 3) + 30yh“H(M, N + 2, 7)] 

+ c[-3H(M.N+ 1,5)+ 15h-H(M, N+ 1,7)]| 

|h(M, N, 3) - 3h%(M, N. 5) + a[3hH(M + 2, N, 5) 

- 1 5h^H(M + 2, N, 7) + 1 2xhH(M + 1 , N, 5) - 30xh^H(M + 1 , N, 7)] 

+ b [3hH(M, N + 2, 5) - 1 5h^H(M, N + 2, 7) + 1 2yhH(M, N + 1 , 5) 

- 30yh^H(M, N + 1 , 7)] + c|9hH(M, N, 5) - ISh^HCM, N, 7)]| 


121 



D.3 CALCULATION OF H INTEGRALS 

In this section we shall compute in closed form each of the integrals 



(D.35) 


P = 


forM= 1,MXQ;N = 1,MXQ-M+ 1 ;K= 1,MXK,2. 

Note that M and N both vary from 1 to MXQ in such a manner that M + N - 1 < MXQ. 

The index K varies from 1 to MXK in steps of 2, which means that K is odd. The 
following values of MXQ and MXK are evident from the previous section. 

Panel Type MXQ MXK 

Source Potential (Flat Panel) 

Source Velocity (Flat Panel) 

Doublet Potential (Flat Panel) 

Doublet Velocity (Flat Panel) 

Source Potential (Curved Panel) 

Source Velocity (Curved Panel) 

Doublet Potential (Curved Panel) 

Doublet Velocity (Curved Panel) 

For the range of indices above, some H(M, N, K) become divergent as the field point 

P(x, y, z) approaches the panel surface S. This is because h approaches zero as P approaches 

S, and therefore p becomes zero when the integration varmbles (5 t?) take on the field 

point values of (x, y). Note that the singularity occurs as P approaches S even though 

the integration variables range over 2. For further analysis see the discussion following 

equations (D.52) and (D.105). 

The integrals H(M, N, K) may be computed with the aid of the following algebraic 
recursion relations. From the definitions of H(M, N, K) and p we can easily derive 
the following identity: 

H(M + 2. N, K) + H(M, N + 2, K) + h"H(M, N, K) = H(M, N, K - 2) (D.37) 
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A second recursion relation can be obtained by considering the identity 

ff ^ |'(£-x)^"2(77-y)’^~ n -wirM - tt i' 

JJ I d^di? = (M - 2) H(M - 2, N, K - 2) 

Z 

-(K-2) H(M,N, K) 

Integrating the left side by parts we obtain: 


4 

(K - 2) H(M, N, K) = (M - 2) H(M - 2, N, K - 2) - ^ F(M - 1, N, K - 2), (D.38) 

I 


and by interchanging the roles of f and 17: 


(K- 2) H(M, N, K) = (N - 2) H(M, N - 2, K - 2) - 


4 

E r'^F(M,N-l.K-2) 


(D.39) 


The summations on the right sides of equations (D.38) and (D.39) are over all four 
sides of 2. 

Here S?is the unit outer normal of the side L (see fig. (D.3)) and F(M, N, K) is the 
line integral for side L, defined by 


de , p= - x)^ + (tj - y)^ + h^ 


(D.40) 


The procedure for evaluating F integrals will be described following equation (D.53). The 
fundamental integrals are H( 1 , 1,3) and F( 1 , 1 , 1 ). Once these two integrals are evaluated, 
the remaining H and F integrals can be evaluated from recursion relations. The details of 
the evaluation of H( 1 , 1,3) are given in section G. 1 of appendix G, and F( 1 , 1 , 1 ) is given 
by equation (D.60). Assuming the F integrals and H(l, 1,3) are known, the recursion 
relations (D.37), (D.38) and (D.39) may be recombined (app. G.2) to yield the efficient 
procedure, given below, for calculating the H integrals. Because some of the H integrals 
are singular on S it is actually necessary to consider three slightly different procedures 
depending upon the relationship of the field point ^ to the panel S. Define djj to be 
the minimum distance of the point (x, y) to the perimeter of 2. If is some small 
number (nominally chosen as 0.01) we have the following three computational procedures. 
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(£ - X, 17 - y) = (ai'j - 


T? 



Figure D.3. - Quadrilateral Geometry 

Procedure 1 : \h\> Sjjdjj, (i.e., P is not “too near” the plane of Z). When the following 

eight steps are performed in the order given, all the H(M, N, K)’s will be obtained for the 
MXQ and MXK given by equation (D.36). 

1 . 4 4 

H(l, 1, l) = -lh|^tan-l [a(C 2 Cj -C 1 C 2 ), ciC 2 + a2fi,£2] +^aF(l,l,l) (D.41) 
1 1 

where the tan’ ^ terms are from H(l, 1,3) (see equation (G.24), and where 

Cj = g- + I hi - + g- , C2 = g^+lhl N/a“ + h“ 

Here tan"^ (y, x) is defined by 0 = tan“^ (sin (t>, cos «^), e[- ir, ir] , and is the same as the 
ATAN2 function of FORTRAN. The sum is over the four sides of S. 
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0 


H(l, 1,K)= 

(K- 2)h^ 


4 

(K-4)H(1, l,K-2)+ X) 5F(1, l,K-2) ; 

1 


K = 3, MXK, 2. 


(D.42) 


3. 


H(2, N, 


1 ) = 


1 

(N+ 1) 


4 4 

h“ j^^Fd, N, 1) +2^aF(2, N, 1) ; 

1 1 


N= 1, MXQ- 1. 


(D.43) 


4. 

H(1,N, 1) = ^ 
N 


4 

-h2(N-2) H(l,N-2, l) + h2^j.^F(l,N- 1, 1) 

1 


4 

+ aF(l,N, 1) 

1 


N = 2, MXQ 


(D.44) 


Note that when N = 2, H(1 , N - 2, 1) need not be computed since it is multiplied by zero. 
The same holds for similar terms in steps 5 and 6. 


5. 


H(M, N, 1) = 


1 

(M + N- 1) 


4 

-h2(M-2)H(M-2, N, D + h^ F N, 1) 

1 


4 

+ 52aF(M,N, 1) 
1 


M = 3, MXQ and N = 1, MXQ - M + 1 


(D.45) 
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6. 


1 


H(l, N, K) = 


(K-2) 


4 

(N-2)H(1,N -2, K-2)- ^ »;^F(1,N- l,K-2) 

1 


N = 2, MXQ and K = 3, MXK, 2 (D.46) 

4 

K-2) : 

1 

N = 1 , MXQ - 1 and K = 3, MXK, 2 (D.47) 


H(2. N,K) = 


1 

(K-2) 


8 . 

H(M, N, K) = - H(M - 2, N + 2, K) - h“H(M - 2, N, K) + H(M - 2, N, K - 2); 


M = 3, MXQ and N = 1 , MXQ - M + 1 and K = 3, MXK, 2 (D.48) 

Procedure 2: I h l< Sj^dn and (x, y, 0) ^ 2, (i.e., ? is “near” the plane of 2 but outside 
the boundary of 2). In this case the recursion defined by equation (D.42) cannot be used 
because h may be zero or close to zero. Hence, H(l, 1, K) for K>3 must be computed 
by some other means. Note that if H( 1 , 1 , K) were known for some lai^e value of K, then 
the recursion of equation (D.42) could be reversed and H(l, 1, K) could be computed for 
successively lower values of K. This is what is done in steps 1 and 2 below. The justifica- 
tion for step 1 is contained in appendix (G.8). 


1. Set 


H(l, 1, NHK + MXK) = 0.0 


where NHK is a positive integer (nominally taken to be 16), and 


(D.49) 


2. Compute j 

H0.1.K-2)= — 

for K = NHK + MXK, 3,-2 

3. All the remaining steps in Procedure 2 (steps 3 through 8) are the same as for 
Procedure 1. 

Procedure 3: |h| < 6j^dj^ and (x, y, 0) e 2, (i.e., ? is “near” the plane of 2 and within 
the boundary of 2). 


r 1 

[h-(K-2)H(l, 1,K)- XlaFd, 1,K-2)J 


(D.50) 
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Define 


H*(M. N. K) = H(M, N, K) - 27ti'(M, N, K) | hl'^+^-K (D 5 ] ) 


where 


/ 0 if M or N is even 
1 for M = 1 


i^(M, N, K)=< 


1 for N = 1 


[1 • 1 • 3 • 5 ... |M - 21] [1 • 1 • 3 • 5 ... |N-2|] 
[(K - 2)(K - 4)(K - 6) ... (K - M - N)] 


(D.52) 


Otherwise. 


and K is odd. 

Then it may be shown (app. G.7) that H* obeys the same recursions as H (except for the 
initial condition (D.41 ) which is irrelevant for procedure 2), and is continuous with respect 
to the field point P when (x, y, 0) e 2. Hence H* may be computed using procedure 2 
(see app. G,8). 

Once H* has been calculated H may be obtained from equation (D.5 1). Care must be 
taken when M + N < K with M and N odd. H* is continuous in the interior of S but 
equation (D.5 1 ) shows that H is singular there when M + N < K with singular part 

(M, N, K) |h|^ + N - K Upon solving (D.5 1) for H and substituting into our previous 
formula for potential and velocity, the only terms involving v that don’t cancel each 
other are of the form I'Th/lhl, where T is continuous. Such terms are responsible for 
jump discontinuities in potential and velocity across S. Recall from appendix C that the 
potential and velocity at a panel can be expressed as upper and lower values or as average 
and difference values (e.g., see equations (C. 5) and (C. 6)). The terms i^Th/|h| make no 
contribution to average values since they are positive on one side of the panel and negative 
with equal magnitude on the other side. Hence, for simplicity we choose to evaluate the 
average value of potential and velocity on S when h = 0 (numerically when |h| < ^Qdpj 
where 6q is nominally chosen as 10" ®). This involves no loss of generality since the 
difference values can be computed directly from equations (C.IO) and (C.15), hence, upper 
and lower surface values can be easily computed from equations (C.9) and (C.13). As a 
practical matter, the desired evaluation will be accomplished simply by replacing h/|h| 
everywhere it appears by signh, where 

( + 1 if h>0 

signh = |-1 if h<0 

'0 if h = 0 (numerically |h| < Sq^H^* 

Note that when h = 0, all terms involving h/|h| will disappear, with the result that average 
values are calculated. 
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We now evaluate the integrals F(M, N, K) for the indices M, N, K required by the H 
evaluation procedures. It is apparent that we need only the following F integrals. 


F(1,1,K) ; 

K= 1,MXFK, 2 




F(M,N, 1) ; 

M= 1,MXQ and 

N= 1, 

MXQ - M + 1 


F(1,N, K) ; 

, N = 2, MXQ and 

K = 3, 

MXK - 2, 2 

(D.53) 


where 


MXFK = 


( MXK-2if Ihl^Sj^d^ 
lNHK + MXK-2if |h| < 


(procedure 1) 

^h^^H (procedure 2) 


(D.54) 


These integrals may be obtained with the aid of three recursion relations. We have the 
following two identities. 

F(M + 2, N, K) + F(M, N + 2, K) + h^FfM, N, K) = F(M, N, K - 2) (D.55) 


and 


i^jF(M + 1 , N, K) + i^^F(M, N + 1 , K) = aF(M, N, K) (D.56) 

Equation (D.55) is the analogue of equation (D.37), and equation (D.56) follows from 
the equation in figure D.2 A third identity is obtained by considering the expression 

E«.,) = 

p 

where p(^, tj) is given by equation (D.40). From figure D.2 we see that ($ - x) and (t? - y) 
are functions only of C along typical side L since a is constant. Thus, we can write 


^ - QE 3(^ - x) ^ 3E 3(t? - y) 
de ” 3({ - x) 0)2 0(r?-y) 
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Performing the indicated operations yields 


-(M- 1,N, K-2) + (N- 1)j^^F(M, N-l.K-2) 

+(K - 2) n^F(M + 1 . N, K) - (K - 2) N + 1 , K) = E(M, N, K - 2) (D.57) 


where 


E(M, N, K) = 


(rj - y)^"^ 




= 4T- 


x)“ + (j? - y)2 + h“ 


(D.58) 


The quantities E(M, N, K) may be evaluated directly or else recursively with the aid 
of the formula 


P(l)=(x2 + xj) P(I-l) -x,X2P(I-2) 


(D.59) 


where 


P(I) = A2X2^-1 - A,x,i-1 


(See app. G.4.) 

The recursion relations (D.55), (D.56) and (D.57) may be recombined to yield the 
efficient procedure for evaluating the required F integrals below. (See app. G.3.) Again, 
the singular behavior of some of the F integrals (near the edges of S) requires a special 
case. Let us define dp to be the minimum distance of the field point ?(x, y, z) to the 
perimeter of S. Then if 6g is some small number (nominally chosen as 0.01) we have the 
following two procedures (see app. G.3). 

Procedure 4: g > 5gdp. (i e., P is not “too near” the perimeter of S). 
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where 


. Note that p = \4^ + 


Here we have used three different (but equal) expressions for the evaluation of F(l, 1, 1) 
in order to avoid possible round off problems from cancellation of negative and positive 
numbers in the argument of the natural logarithm. 


F(l, 1,K) = ^5 ^ — [(K-3)F(1, 1, K-2)-»'„ E(2, l.K-2)] 

5 Ed, 2, K- 2)1; 



K = 3, MXFK, 2. (MXFK is defined in equation (D.54). (D.6 1 ) 

3. a. If |v^| < 

(i) F(1,N, 1) = ,-;^, r(2N-3)a«/_F(l,N-i, 1) 

(N - 1) L ' 

- (N - 2)(a2 + 12 ^ 2 ^ 2 ) f( 1,N-2, 1) + vjE(l,N- 1,- 1)] ; 


N = 2, MXQ 


(D.62) 
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M = 2, MXQ and N = 1 , MXQ - M + 1 


(D.63) 


b. If < |.,| 

(i) F(M, 1, = [(2M-3)S^'^F(M- 1, 1, 1) 

+ v^\'^) F(M-2, 1, 1,- 1)] 

M = 2, MXQ 


(ii) 


- Vt 


5 a 

F(M, N, 1) = — ^F(M+ 1,N- 1, 1) + — F(M, N- 1, 1); 


(D.64) 


N = 2, MXQ and M=1,MXQ-N+1 


(D.65) 



4 . 


E(l, l.K-2) 


F(l,2, K) = »'^aF(l, 

K = 3, MXK-2, 2 


(D.66) 


5- F(1,N, K) = 2ai'^F(l,N- 1,K) - (a^ + F(l, N - 2, K) 

+ i^^“F(l,N-2, K-2) : 

N = 3.MXQ and K = 3, MXK-2, 2 


(D.67) 


Procedure 5: g < 6gdp, (i.e., P is close to the perimeter of S). 

, (See app. G.8) (D.68) 

F(l, 1,MXFK + NFK) = 0 

where NFK is a positive integer (nominally taken as 16) 

F(l, 1 , K - 2) = ^ rg“(K-2) F(l, l,K) + i'^E(2, 1, K - 2)-p^E(1, 2, K- 2)] ; 

(K • 3) L 

K = MXFK + NFK, 5, 2 (D.69) 

3. F ( 1 , 1 , I ) as well as other F integrals may be computed in the same manner as for 
procedure 4. 

This completes the calculation of the H integrals. 
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D.4 EVALUATION OF SOURCE AND DOUBLET INTEGRALS 
FOR A DISTANT FIELD POINT 

If the field point P is a large distance from S the approximation (D.I8) may be replaced by 
an approximation based on this fact. Let 

P = |P| and Q = |Q| (D.70) 


Then 


Let 


Then 


]_ 

R 


I 



- 2(P • q)+ q2 


T 


1 

I ^ - 2(P • q)^ 

J p2 


(D.7I) 


1 Max 

P (?,Tj)eS 


{Q} 


(D.72) 


- 2(p 


q)+q2 


< 2e + e" 


(D.73) 


Hence, if 


e < < 1 


(D.74) 


we have 


1 


1 , k(p . q) K 1 
pK pK + 2 2 pK + 2 




(D.75) 
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Only the first three terms of the binominal expansion are displayed (monopole, dipole, and 
quadrapole). In practice this expansion is used only when e is less than 1/5. All three 
terms are used unless e is less than 1/8 in which case only the first two terms are required. 

Substituting (D.75) into (D.IO) and using hypothesis (D.6) we obtain for the source potential 


^ = OqI( 1, l) + Ofl(2, l) + a„I(l,2) 




(D.76) 


Here 


KM, 



1 

3 /a - *\ I 

(E2 - p)+ - 


E3 +-(P ■ E4P) } 



2' L\ 


(D.77) 


where 

El = C(M, N) 


(D.78) 


E2 = [C(M + 1 , N), C(M, N + 1 ), aC(M + 2, N) + bC(M, N +2)] 


(D.79) 


E3 = C(M + 2, N) + C(M, N + 2) 


(D.80) 


E4= /C(M + 2, N) 

C(M + 1 , N + 1 ) 
aC(M + 3, N) + 
bC(M + 1, N + 2) 


C(M+ 1,N+ 1) aC(M + 3, N) + bC(M + l,N + 2)' 

C(M, N + 2) aC(M + 2, N+ l) + bC(M, N + 3) 

aC(M + 2, N+D+ 0 

bC(M, N+3) 


KD.81 


A 

P 


-► 

P 

P 


(D.82) 
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and 


(D.83) 


C(M, N) = 



,N-1 


d{ dTj 


In computing the elements of the matrix E4 hypothesis (D.6) has been employed. The 
integrals C(M, N) will be computed later (see equation (D.94)). The induced source 
velocity may be obtained by differentiating equation (D.76), and is given by 


v = OqJ( 1, l) + a^J(2, l) + o,^J(l,2) 

where 



(D.84) 


(D.85) 


A similar expansion may be obtained for the doublet induced potential and velocity. 
Substituting the approximation (D.75) into equation (D.3) and using hypothesis (D.6) 
we obtain, 


d> = H0 1(1, l) + p^ 1(2, D + p^ I(T,2)+^Pj^ 1(3, 1) + Mj^I(2, 2) 1(1, 3) 


Here 


I(M,N) = — 
4ir 





(D.86) 


(D.87) 


where 


E2 = [- 2aC(M + I , N), - 2bC(M, N + 1 ), C(M, N)] 


(D.88) 


E3 = - aC(M + 2, N) - bC(M, N +2) 


(D.89) 
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and 



-(a + b)C(M+ 1,N + 1) 


(a + b)C(M+ 1.N+ 1) -2bC(M, N + 2) 


1 

-C(M, N + 1) 


1 

;C(M + 1,N) 

^C(M,N+ 1) 

aC(M + 2, N) 
+ bC(M. N + 2) 


(D.90) 


Note that these expressions differ from the E2, E3 and E4 given above for the source. 

Again, hypothesis (D.6) has been used in the computation of the matrix E4. The induced 
velocity may be obtained by differentiating equation (D.86) and is given by 



(D.91) 


Here 


J(M, N) = — 
47T 





r 1 


1 

— ► / — ► a\a 

1 



E2-3 E2 • P P 

+ 

3E3 P - 1 5( P • E4P)P + 6E4P 


p3 

\ ' 


V / 



(D.92) 


Note that the quadrapole term of expansion (D.75) is not used for doublet panel induced 
potential and velocity. This is primarily because of the complexity of this term and the 
resultant fact that evaluation of this term is only marginally more efficient than evaluation 
of doublet velocity and potential from the formulas of the previous sections. 

The computation of the C(M, N) integrals of equation (D.78) follows from the computation 
of the H integrals of section (D.3) by noting that C(M, N) = H(M, N, 0) with x = y = 0. 

The range for the indices M and N is the same as that for the H integrals, i.e., 

M=1,MXQ; N=1,MXQ-M+1 (D.93) 

where MXQ is given in (D.36). Setting K = 2 and x = y = 0 in equations (G.42) and (G.43) 
and adding these equations, we obtain 

4 

(M + N) C(M, N) = X + 1, N, 0) + F(M, N + 1, 0)] 

1 

Upon substituting equation (D.56) on the right we obtain 
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C(M, N) = — — — XI aC(M, N); M = 1 , MXQ and N = 1 , MXQ - M + 1 (D.94) 

(M + N) , 


where 


G(M, N) = F(M, N, 0) (evaluated at x = y = 0) 


= J fM-' ,N-I 


dC; M=1,MXQ and N = 1, MXQ - M + 1 (D.95) 


The integrals G(M, N) may be obtained with the aid of two recursion relations. By setting 
K = 0 and x = y = 0 in equation (D.56) we have 

v^G(M + 1 , N) + J^^G(M, N + 1 ) = aG(M, N) (D.96) 


By setting K = 2 and x = y = 0 in equation (D.57) we obtain 

-(M - 1 ) v^G{M - 1 , N) + (N - 1 ) J^jG(M, N - 1 ) = D(M, N) (D.97) 


where 

D(M, N) = E(M, N, 0) (evaluated at x = y = 0) 

2 (D.98) 

= jM - 1 - 1 

1 

The quantities D(M, N) may be evaluated directly or else recursively with the aid of 
equation (D.59) (see app. G.5). The recursion relations (D.96) and (D.97) may be recombined 
to yield an efficient procedure for evaluating the G integrals, (see app. G.5). 

Procedure 6; 


1. a. 

If 

(i)G(l,N) = ^D(l,N+l): N=1,MXQ 

« (D.99) 
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(ii) 

% a 

G(M, N) = -— G(M - 1,N + 1)+ — G(M- 1,N) ; 
M = 2, MXQ and N = 1 , MXQ - M + 1 


(D.lOO) 


b. 


If < \Pj_ 


(i) 1 

G(M, 1) = D(M+ 1, 1); M = 1,MXQ 

Mu 


(D.IOI) 


(ii) a 

G(M, N) = -— G(M + 1, N- 1) + — G(M, N- 1); 

V 

‘^rt V 


N = 2, MXQ and M = 1 , MXQ - N + 1 


(D.102) 


This completes the evaluation of the source and doublet potential and velocity for a 
distant field point. 

D.5 BEHAVIOR OF INDUCED POTENTIAL AND VELOCITIES 

In this section we shall study the behavior of the potential and velocity induced by source 
and doublet panels. We shall restrict our attention to flat panels since the addition of 
curvature produces no qualitative change in behavior. First, consider the potential ^ 
induced by a flat source panel. Setting a and b zero in equations (D.IO), {D.15) and (D.16) 
we obtain 


<P = 



(D.103) 


where 

R= Jii - x)2 + (tj - y)2 + h^ =p 


and 


h = z 
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We note from (D.103) that 0 is regular everywhere in finite space except perhaps on Z. 
From equation (D.75) we see that near infinity 


1 





Z 


(D.104) 


Since 0 vanishes as 1/P near P = «>, 0 is regular at infinity as well. To examine the behavior 
of 0 near Z we use equations (D.2C^ and (D.22) with a = b = 0. 


Then 


0 = o(x,y)[-^H(l, 1, 1)] 

+ Ojj(x, y) H(2, 1, l)J + ay(x, y) H(l, 2, 1)J (D.105) 

Appendix G.6 shows that h"^ H(M, N, K) is continuous near Z and equal to zero on Z if 
J + M + N > K, and h'^ H(M, N, K) is bounded there if J + M + N = K. Consequently, 0 is 
continuous everywhere, in particular on Z. 

From equation (D.103) we have 


V0 = 



dfdr? 


(D.106) 


Equation (D. 106) shows v is regular every where in finite space except perhaps on Z. 
From equation (D.104), we see that near infinity 




so that"v = 0(1 /P^) if total source strength is zero and'^= 0(l/P^) otherwise. 


(D.107) 


The behavior of v near Z will now be examined. For this purpose one can use equations 
(D.27) and (D.28) with a = b = 0. However, it is more instructive to derive an alternate 
expression forV! Differentiating equation (D.103) with respect to x we obtain 


3x 




(D.108) 


139 



Integrating by parts we obtain 


3x 




dC 


(D.I09) 


Similarly 




Finally 


30 

3z 


/•(=) 


d^dTJ 


Combining all three equations we have 

30 30 


_/30 30 30 

\ 3x ’ 3y ’ 3z / ^ 




where 




2 


(D.llO) 


(D.lll) 


(D.l 12) 


(D.113) 


(D.l 14) 


(D.l 15) 


(Note that here is not the same quantity as the average velocity of equation (C.5).) 
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First consider the vector From equations (D. I 1 3), (D. 8 ), (D.22), (D.25), and (D. 1 5) 
we have 




(x, y), ay(x, y). 



H(U 1, I) 


(D.l 16) 


Again H(l, 1, 1) is continuous everywhere, hence v^ is continuous everywhere. 

Consider Vg for a typical side L of S. Note first that if L is a common edge of X and an 
adjacent panel across which source strength is continuous and the surface slopes are continuous 
(i.e., 0) of L = - 0) of the adjacent panel) then the behavior of vg is irrelevant 

since "vg is cancelled by the same component of the adjacent panel induced velocity. In any 
event, the behavior of\Tg can be established as follows. From equations (D.l 14), (D.22), 

(D. 15), and (D.40) we have 




, ,F(1, 1,1), ^ ^F(2, 1. 1) . ^ ,F(1,2, 1) 

a(x, y) ^ + a^(x, y) 7 - + a,,(x, y) 


47T 


47T 


Att 


(D.l 17) 


In appendix G .6 we show that g*^ F(M, N, K) is continuous everywhere ifJ + M + N>K + 1. 
Therefore, any discontinuity in 7g must arise from F(I, 1, 1). Again, from equation (D.40) 
F( 1 , 1 , 1 ) is continuous everywhere except possibly on L, where g = 0. From equation 
(D.60) we derive the following asymptotic formulas for F(l, 1 , 1) as F approaches a point 


Pl6L. 


F(l, 1, 1) 


ln(l/g^) ; in the interior of L (at £) 
ln(l/p); C| • ^2 ^ 

^ - w 

V ln(p/g”) ; ^2 ■ ®1 and at endpoint of L 


(D.l 18) 


Here p = |?- Pj^l = s/g^ + . We see that has a logarithmic singularity on L 

which is proportional to local source strength. 


Finally, we obtain from equations (D.l 15), (D.22) and (D.25) that 


vc= (0,0, 1) 


a(x, y) 


hH(l, 1, 3) 


+ a (x, y) 

^ An 


47T 

hH(2, 1, 3) 


. , ,hH(l,2,3) 

+ Oy(x,y)^^^ . 


(D.l 19) 
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Any possible discontinuity in must arise from the term ^ can only occur 

on S. From equation (G. 24) we have 

hH(l, [a(62C, -6,C2), c,C2 + a2fijC2] (D.120) 

4ir 47r\|h|/ | 

One can show (see sec. G.l of app. G) that the sum of arctangents on the right is continuous 
everywhere except on the perimeter of Z where it remains bounded. Moreover, when 
h 0 this sum approaches the value lit inside Z and 0 outside Z. Consequently, as h ^ 0, 

approaches 0 outside Z, 1/2 on the upper surface of Z and -1/2 on the lower 
surface of Z. ^ ^ defined to be zero on Z (actually when (x, y) e Z and |h| < 

where 6q is nominally chosen as 10'^) and this results in the computation of an average 
velocity on Z as discussed in section D.3. The discontinuities of v^ are essentially the same 
except for the proportionality factor a(x, y). Note that this behavior gives the tirst term on 
the righthand side of equation (C.IO). 

The behavior of ^ in finite space can be summarized as follows. The normal component 
of 7 is bounded every where but discontinuous on Z. As the field point approaches the 
source panel plane, this component is zero outside Z, equal to 1/2 local source strength 
on the upper side of Z and -1/2 local source strength on the lower side of Z. (This component 
is defined to be zero on Z.) The tangential components of v are continuous every where 
except on each edge L of Z. As the field point F approaches a point interior of 

L, T has the following singular behavior due to the characteristics of vg: 

-2alog(g)^®n (D.121) 

~ 4ir 

a — ^ a 

Here o is the source strength at P^, n is the normal to Z at P^, C is a unit vector along 
L such that points out of Z and g is the distance from ^ to the line containing L. 
If ^ approaches a point at a comer of Z, the singular behavior of 7 is derived by 
summing the singular contributions of the two intersecting edges. These contributions are 
described by the right side of equation (D.121) with but slight modification to the factor 
-2 log(g) arising from the alternate expressions of equation (D.l 18). 

Next we consider the potential 0 induced by a flat doublet panel. Setting a and b zero 
in equation (D.3) we have 

0 

Z 

We note from (D.l 22) that 0 is regular everywhere in finite space except perhaps on Z. 
From equation (D.75) we see that near infinity 




\4wR^/ 


(D.l 22) 
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(D.123) 


2 

so that 0 is regular at infinity as well. 

To examine the behavior of 0 near 2 we use equations (D.33) and (D.34) with a = b = 0. 
Then 


r h 1 

r h ^ 1 

M(x,y) — H(l, 1,3) + 

M,(x,y)^-H(2,l,3)J 

fh 1 

1 f h 1 


+ H<3,1,3)J 

r h 


+ H^y(x,y)]^H(2, 2,3) 

+ -Myy(x,y)|^- H(l, 3,3)j 


(D.124) 


Any possible discontinuity in 0 must arise from the first term on the right. We have 
already analyzed connection with equation (D.120). Consequently, we 

47T 

can say that 0 is bounded in a neighborhood of 2. Moreover, as the field point 
approaches the doublet panel plane 0 is zero outside 2 equal to 1 /2 local doublet 
strength on the upper surface and -1/2 local doublet strength on the lower surface of 2. 
(0 is defined as zero on 2.) 


From equation (D.122) we have 


v = = 


=/T4— 

JJ L 4jrR- 


O ^ 3hR 


4jtR^ 


d^dT? 


(D.125) 


Equation (D. 1 25) shows that v is regular everywhere in finite space except perhaps on S. 
From equation (D.123) we see that near infinity 


7; 


(0, 0, 1) 3hP 
d3 d5 




(D.126) 


so that "v = 0(1 /P^) if total doublet strength is zero and"v = 0(1 /P^) otherwise. 


The behavior of “v near 2 will now be examined. Differentiating equation (D.122) with 
respect to x we have 
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Similarly 


Finally 



30 

3y 




(J-x) ^ b^l (rj-y) 


dfdTj 


JJ 47tR^ 47tR^ 

2 

L £ 

= -/*^^de- /■fe, J-L d£ 

y 47 tR^ y ^ 4 itR 


im 


brp- 


47tR 


dfdt7 


(D.127) 


(D.128) 


dfdr? 


(D.129) 
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Combining all three equations we have 


V = 


^ ^ ^ 

3x ’ 9y ' 9z , 


= ^A+ E Vg+ E Vc + Vj) 


1 


1 


(D.130) 


where 


// 

2 

/ 


3^/l/ 3^ 


ar ar? 


arj^J 


g = # (- hi^^, - hi'Tj- ■ a) M 


47tR- 


47tR 


dC 


dfdTj 


-■//(S'S ")i - 


(D.131) 


(D.132) 


(D.133) 


(D.134) 


(Note that vgj is not the same as the difference velocity of equation (C.6).) 

First consider the vector v^- From equations (D. 131), (D.9), (D.25), and (D.3 1 ) we have 


= (0, 0, l)|^Mxx(’‘-y) + Myy(x, y)j 1^^ H(l, 1, l)j 


(D.135) 


Again H(l, 1, 1) is continuous everywhere, hence v^^ is continuous everywhere. 

Next consider vg for a typical side L of E. Note first that if L is a common edge of E and 
an adjacent panel across which doublet strength is continuous then the behavior of is 
irrelevant since vg is cancelled by the same component of the adjacent panel induced 
velocity. In any event the behavior of Vg can be established as follows. From its definition 
vg is clearly continuous except possibly on L. To examine the behavior of 7g on L we 
have from equations (D. 132), (D.9), (D.3 1 ), and (D.40) that 
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gF(l, 1,3) 
■y) 4n 


^ ^ gF(2, 1,3) , ^ gF(l,2,3) 

+ +My(x,y) 

1 gF(3,l,3)^ ^ ' gF(2, 2, 3) . 1 , 

+T/^vv(x, y) ::7 +Mxv(’‘’y) 4_ 


(D.136) 


The direction of vr is the unit vector 


'Jm 

g ’ g ’ g / le ® (?L 


(D.137) 


A “► — 

Here £ is the unit vector Ut, 0) along L and Pl is the closest point on L to P. The 
magnitude of is the term in the square brackets on the right side of equation (D.136). 
The second and third terms are bounded and the last three terms are continuous in a neigh- 
borhood of L. However, the first term is unbounded as approaches L. From equation 
(D.61) we have the following asymptotic formulas for the coefficient of p(x, y) as P 

approaches 


gF(l,l,3) ~ ; Pt in the interior of L 

g ^ 


. 9.2 >0 and Pj^ at endpoint of L 


1 + J\ - g^/p“ 


1 + n/i "g^/P“ : Cj • £9 < 0 and at endpoint of L 


(D.138) 


Consequently vg becomes unbounded as the reciprocal of the distance from P to L. 

Next we consider v’^- for a typical side L of 2. Note first that if L is a common edge of 2 
and an adjacent panel across which the derivative of n perpendicular to L and surface slope 
are continuous then the behavior of "v^ is irrelevant since v^ is cancelled by the same 
component of the adjacent panel induced velocity, "v^; is clearly continuous except perhaps 
on L. From equations (D.133), (D.9), (D.31), and (D.40) we have 
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VC = (0, 0, - 1 ) j [mx(x, y) + /iy(x, y) I' J — — — 


+ [Mxx(x- y) «^J + Mxy(x, y) I'ri] 


F(2. 1, 1 ) 
47T 




xy(x, y) l'j+Myy(X,y) t',,] 


F(l,2, 1) 1 
47T 


(D.139) 


Any discontinuity in must come from the factor F(l, 1, 1). The asymptotic behavior 
of this function is displayed in equation (D.l 18). We see then that the normal component 
of 'vq has a logarithmic singularity on L proportional to the local derivative of m in a direc- 
tion perpendicular to L, 

Finally, we consider Clearly is continuous except perhaps on 2. From equations 
(D.l 34), (D.9), (D.31), and (D.24), we have 





+[Mxx(x,y), 

+ [Mxy(x, y) , 


Mxy(x, y). 

Myy(X, y), 


hH(2, 1,3) 


•][ 

.j[-^ 


3) 


] 

] 


(D.l 40) 


The quantities hH(2, 1, 3) and hH(l, 2, 3) are continuous everywhere, however, 

_ hH(l, 1, 3) is discontinuous on S as discussed previously (equation (D.l 20)). The behavior 
^■Tr ^ 

of is essentially the same except for a proportionality factor vm(x, y). 


The behavior of v in finite space can be summarized as follows, v is continuous everywhere 
in finite space except perhaps on 2. Asa field point ? approaches a point the interior 
of 2 the normal component of "v is continuous at P^but the t^gendal components 
approach l/2v/i at?2 ^ approaches 2 from above and -1/2^M at if P approaches 2 
from below. (The tangential components of ^ are defined to be 0 on 2 so that they are 
zero everywhere in the plane of 2.) Note that it is this behavior that gives the second term 
on the righthand side of equation (C.IO). 


As a field point P approaches a point Pl in the interior of an edge L of 2, v has the 
following singular behavior due to the characteristics of and v^^ discussed above: 
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(D.141) 


2n 


~ n 

£®g 

4jrg 


2C 


(n « V ji) 

47T 


log(g) n + bounded terms 


Here /i is doublet strength and vm doublet gradient at P^. n is the unit normal to Z ^ , 
8 is a unit vector along L such that 8® n points out of S and g = gg is^he vector from r 
to the closest point on the line containing L. If ^ approaches a point at a corner of Z, 
the asymptotic behavior of v is derived by summing the contributions of the two inter- 
secting edges. These contributions are described by the right side of equation (D.141) with 
but slight modification to the factors 2/g and 2 log(g), arising from the alternate forms in 
equations (D. 138) and (D.l 18) respectively. 

D.6 DERIVATION OF BOUNDARY VALUE PROBLEM INFLUENCE 
COEFFICIENT EQUATIONS 


In section D.2 we derived the expressions for perturbation potential and velocity due to 
source or doublet panels. These expressions are in terms of the panel singularity distribution 
coefficients of equations (D.8) and (D.9). Before boundary conditions can be imposed, _ 
these expressions must be re-expressed in terms of the unknown singularity parameters X of 
appendix B. In this section we describe how this is done for the specific case of the 
velocities due to a doublet/analysis network. 

From equations (D.33) and (D.3 1 ) we obtain 


where 


V = 


LTJ 




M 

TflJ 


(D.l 42) 


LU = 


j(i, 1) 


xj(l, 1) 
+ J(2, 1) 


yJ(i, 1) 
+ J(1,2) 


-x2j(i, 1 ) j xyj(l, 1) 

I ^ 

+ xJ(2, 1) I +yJ{2, 1) 

1 I -► 

+ -J(3, 1) I +xJ(l,2) 


I 

I +J(2, 2) 


1)1 

+ yJ(i,2) 

1 -► 

+-J(1,3) 
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Using equation (D.34), the scalar form of equation (D.142) is 


"x 


1) I 

X J^d, 1) + Jx(2, 1) 

1 

i 

\ 

— ► 

• = 

Jyd, 1) 1 

X Jy d, 1) + Jy(2, 1) 

1 

1 etc. 

1 

— ► 

''z 


Jz(M)| 

X J^(l, \) + }^(2, 1) 

1 

1 

1 



^ 1 


1 M 


n-n 

^1717 


(D.143) 


Equation (D.142) gives the velocity at field point i (i.e., at Xj, Vj, Zj) due to the doublet 
distribution 


1 ^ 1 o 

ri) = HQ + + + MJt? f^7+ — 7*7777 

associated with a particular panel. If we label this panel as panel number k, we rewrite 
equation (D. 142) as 




(D.144) 


From appendix B, recall that the six unknown doublet strength coefficients mq 

panel k are expressed in terms of a subset of the unknown doublet singularity parameters 

X; denote this subset as Xj. From equations (B.3) and (B.6) we have 



(D.145) 


where Nj^ is the number of doublet singularity parameters associated with panel k. 
Equations (D. 1 44) and (D. 1 45) give 


%(j)^= Wk Mk Mk 


(D.146) 
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where the subscript (j) is used to indicate that the velocity due to the single panel k depends 
on several singularity parameters Xj. 

Now imagine that the network contains nine panels (e.g., as in the type 2 schematic of 
figure Then the velocity at point i due to these nine panels, and the free-stream 
velocity is 


9 

""1 = ^00+ = Vj (X) (D.147) 

k= 1 

where Vj (X) denotes that the velocity is now in terms of all the singularity parameters 
associated with the network. For the case of figure B.2, type 2, there are 25 singularity 
parameters. 

Next, let the field point i be one of the 25 control points on the network (see fig. C.l, 
type 2). Imposing the boundary condition = 0 (where n^ is the unit normal at control 

point i) at each of these control points gives 

25 

E V'"i ' ■'^-•"1 

k= 1 

When cast in matrix form, equation (D.148) becomes 


25 X 25 

h] 



(D.149) 


Each row i of the aerodynamic influence coefficient matrix [Ay] represents a boundary 
condition imposed at one of the 25 control points. Each column j corresponds to one 
of the 25 singularity parameters in |x|. The matrix [Ay] is constructed one row (control 
point) at a time. For each row, one cycles through the panels and enters the contributions 
of each panel to the appropriate columns of [Ay] . 

For more than a single network, the procedure is exactly the same except that the matrices 
in equation (D.149) expand in size so as to incorporate all the panels, all the singularity 
parameters, and all the control points of every network. The general form of equation (D.149) 
is then 


M x M M X 1 M X 1 
[A] {X[ = jb( 


(D.150) 
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where M is the toUl number of singularity parameters (and control points) for all the 
networks. Hence |x| can be solved for, and tlien the resulting velocities can be calculated 
from equation (D. 147). (The value k = 9 appearing in equation (D.147) would be replaced 
with the total number of panels in all the networks.) Knowing the velocities, the pressures, 
forces and moments are then computed as shown in appendix F. 
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APPENDIX E 
EQUATION SOLUTION 


In this appendix we give a brief description of the numerical method used in the pilot code 
to solve the influence coefficient equation set. This equation set is denoted in matrix form 
by 

AX = B (E.l) 

where A is an NxN matrix of coefficients, B is an NxM matrix of rightjiand sides and X 
is an NxM matrix of unknowns. The method assumes that the matrix A resides in 
mass storage in such a manner that ea^ row is a record which can be accessed in a 
sequential manner only. The matrix B is also assumed to be stored in* the same way and 
once the matrix X has been computed it is likewise stored in the same way. 

For large values of N, the data transfers to and from mass storage during solution became 
excessive if the row storage format is retained. Hence, the metho^begins^y reformatting 
the A and B matrices in block form. Specifically, the matrices A and B are partitioned 
into blocks (e.g., fig. E.l) 


>1 

II 

>1 

i= l,n 

and j = 1 , n 

(E.2) 

B = Bij ; 

i = 1 , n 

and j = l,m 

(E.3) 


which are stored as randomly accessible records. The values of n and m are determined 
by the available central memory. 

The Grout decomposition algorithm is then used to solve equation (E.l). This algorithm 
employs the substitution 

lu = A (E-4) 

to obtain 

IY = B (E.5) 

where 

UX = "V (E.6) 

The matrix L is lower triangular and U is upper triangular normalized by ones along the 
diagonal. The decomposition of equation (E.4) is accomplished as follows. The matrices 
L and U must satisfy the relationship 

N 

Eeik“kj="ij (E.7) 

k= 1 
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where Cjj, uy, ay are the elements in the ith row and jth jolumn of L, iJ and A respectively. 
From this equation we obtain 


and 


j - 1 

^ ®ik^kj ’ 

k= 1 


i - 1 

''ij^ (^ij‘ S ) Ai ' ' ^ 

k= 1 


(E.8) 


(E.9) 


Since A is physically stored as its component submatrices Ay, the operations of equations 
(E.8) and (E.9) are actually performed by forming the submatrices 


min(i, j) - 1 


(E.IO) 


k= 1 


and then solving 


and decomposing 


Lij^jj 

= ^ij 

for 

^j 

when 

i>j. 

LiiUij 

= CiJ 

for 

UiJ 

when 

i < j 



Lii Uii 

when 

i = j 



(E.ll) 

(E.12) 

(E.13) 


The forward substitution of equation (E.5) is accomplished as follows. The matrix Y must 
satisfy the relationship 


so that 


^ *^ij ^kj *^ij 

k= 1 



(E.14) 


(E.15) 


where yy, by are the elements in the ith row and jth column of Y and B, respectively. 
The block operation is performed by first forming 
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(E.16) 


i- 1 _ 

Dij - Sy- S 

k= 1 

and then solving 


Lii Yij = Dij 


(E.17) 


by forward substitution. 

The backward substitution of equation (E.6) is accomplished as follows. The solution X 
must satisfy the relationship 


n 

2] “ik^kj = yij 
k= 1 

where Xy is the element in the ith row and jth column of X. Since 


we have 

n 

^ij “ yy ’ E ^ik ’‘kj 

k = i+1 

The corresponding block operation is performed by first forming 

Eij = Yy - f: UikXkj 
k = i+l 


(E.18) 


(E.19) 


(E.20) 


(E.21) 


and then solving 


^ii ^ij ^ij 


(E.22) 


by backward substitution. 

The present method does in-block partial pivoting. The pivoting is established during the 
in-core decomposition step. That is the quantity 


_ _ * ■ 1 _ _ 

^ii “ ^ii " ^ Mk ^ki 
k = 1 


(E.23) 
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if formed without any pivoting considerations. The decomposition of Cjj is performed with 
pivoting, that is 


C = P L • U (E-24) 
^11 ^1 hi 

where Pj is the identity matrix with row interchanges corresponding to the required pivoting. 
At this point, the identical row swaps in Ly for j < i and By for j = 1, m could be performed. 
However, this is not done. Instead, Pj is stored as a vector describing the row swaps required 
to generate Pj from the identity matrix. Each time Ljjis used to perform forward substitu- 
tions, the swaps dictated by this vector are first performed on the righthand side. 


The method used to detect singularities in the system uses discontinuity of the pivot elements 
to flag the singularity. Each proposed pivot element is compared in absolute value to the 
minimum previously accepted pivot element. If the proposed element is much less in 
magnitude, a singularity is declared. While this method is adequate for the detection of 
obvious singularities (row identically zero, two rows the same, etc.), it does not adequately 
detect ill-conditioned problems. An additional check is provided to detect ill-conditioned 
systems. The method approximates error growth by the following scheme. Let 


0 


(E.25) 


and 


^ij, k ^ij, k-1 * ^ik, k-1 ” ^ik, k-1 ^kj, k-l/^kk, k“l (E.26) 

The relative error e is approximated by 


X io->4 (E.27) 


The quotient ay ratio of the remaining matrix elements at elimination step 

k divided by the kth piVot. The product of this ratio with the growth term (k + 1 ) k gives 
the error growth ratio. The assumed initial error is 10”^^. A singularity is flagged when e 
exceeds 10”^. 


Max 

e= 1 <k<N 
i,j>k 


k(k+ 1) 


Pij.kl 
l^kk, kl 
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APPENDIX F 

COMPUTATION OF AERODYNAMIC QUANTITIES 

In this appendix we describe the manner in which certain local and global quantities of 
aerodynamic interest ^e computed. Our point of departure is the computation of the 
average total velocity at each panel center control point along with the distribution of 
singulariW strength on each panel. From the singularity distribution the total difference 
velocity Vp may be calculated from equation (C.IO) and then the total upper surface velocity 
Vjj and total lower surface velocity computed from equation (C.9). (However, 

note the discussion following equation (C.IO) which implies that in the case of a source 
panel superimposed on a doublet panel the normal difference velocity is computed correctly 
only at the source panel control point, whereas the tangential difference velocity is computed 
correctly only at the doublet panel control point.) 

The two velocities Vy and are of most interest in aerodynamic applications and all 
aerodynamic quantities are computed separately for each type of velocity. Hence, let V 
be one of these velocities. We then define a pressure coefficient Cp at each panel center 
control point by the formula 

Cp= 1 -(v . v)/(v^. vj (F.l) 

where is the freestream velocity. (Because of limitations in the existing pilot code logic 
the computed Cp will be in error for control points on superimposed source and doublet 
sheets, however, in the usual case of small or zero normal velocities the value of Cp at the 
doublet control points will be quite accurate). 

A linear distribution of Cp on each panel can be calculated from the values of Cp at the 
control points. The method used is identical to that for computing the linear source distribu- 
tion of network type number 1 (see app. B). On any panel S, we then have 

Cp (f, V) = Cp^ + Cp^ ^ + Cp^ V (F-2) 

The force coefficient vector Cp on S is defined by 

Cp = ^ Jjf Cpnds (F.3) 

^ S 

where Sj^ is a specified reference area. 

Substituting (F.2) into (F.3) and using (D.30) and (D.12) we obtain the following 
expression for ?p in local coordinates, 
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Cp = ^ JJ (Cp^ + Cp^ f + Cp^ T?) (- 2a{, - 2btj, I ) d{di7 

= ^ Cp [-2aC(2, l),-2bC(l,2),C(l, 1)1 

Sr PoL J 

+ ^ Cp [- 2aC(3, 1 ), - 2bC(2, 2), C(2, 1 )]+ ^ C_ [- 2aC(2, 2), - 2bC( 1 , 3), C( 1 , 2)1 

Sr e Sr t'T? 


(F.4) 


where C(M, N) is defined by equation (D.83) and computed in appendix IX section 4. 

(Cp may be expressed in global coordinates via a premultiplication by A ^ where A is the 
transformation matrix defined in equation (A. 17).) The sum of the force coefficient vectors 
for each panel in a network yields a force coefficient for the whole network. The force 
coefficient vectors for any collection of networks can be combined to yield a force coefficient 
vector for the configuration represented by that collection. This force coeffident vector 
can then be resolved into components along (drag) and perpendicular to (lift and 
side force). 


The moment coefficient Cj^ at a point Rj^ about an axis t|^ is defined by 

1 






CpdS 


(F.5) 


where Tr is a reference length and R is the integration point on S. Equation (F.5) is 
equivalent to 


C 


1 A 

M”^^R ‘ 
Tr 




(F.6) 


where is the center of the local coordinate system on S. 


Here Cp is defined by 


Ct: - 


^R 


II 


Q ® n Cp dS 


(F.7) 


where 

Q=R-Ro 

^ — ► 

Cp may be computed in local coordiantes in the same manner as Cp, i.e., 
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2 

" % If (^*’° ^ ^ ^ ’^)] 

2 

= ^ So[^^‘’ 2)’-C(2, 1) , 2(a-b) C(2, 2)] 


Cp [C(2, 2), -C(3, 1), 2(a-b)C(3,2)1 


+ 7^ Cp [C(l,3), -C(2, 2), 2(a-b)C(2,3)1 

Sr Pt/L J 


(F.8) 


Here we have neglected terms of second order in a and b on account of hypothesis (D.6). 
The vector ?£ may be expressed in global coordinates in the same way as ?p. 
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APPENDIX G 

DERIVATION OF RESULTS GIVEN WITHOUT PROOF 
IN APPENDIX D 

In Appendix D, several results were given without proof. This was done in order to emphasize 
the procedures required to generate the influence coefficients. In this appendix we provide 
the proofs and derivations that are missing from Appendix D. 

G.I. EVALUATION OF H(l, 1,3) 

In section (D.3), H( 1 , 1,3) was referred to as the fundamental H integral, since all the other 
H (M, N, K) integrals are obtained from it and F(l, 1, 1) by the recursions given in 
procedures 1, 2, 3, 4, and 5. In this section, we give a detailed derivation of the closed form 
integration for H(l, 1, 3), and then describe the behavior of hH( 1 , 1, 3). 

From equations (D.25), (D.15) and (D. 16), we have 


where 


H(l, 1,3) = 


J 


d^di7 

„3 


(G.I) 


p = + h^ 

r= s/(t - x)^ + (t? - y)2 

h = z-ZQ (G.2) 


are illustrated in figure D.2. In equation (G. 1 ), the integration is over the flat surface 
dS = dfdrj and h is a constant as far as the integration is concerned. Changing to polar 
coordinates, equation (G. 1 ) becomes 


H(l, 1, 3) = 



(G.3) 


where the upper limit on r extends to the boundary of 2. The geometry of the situation 
is shown in figure G.I, where (x, y,0 ) is the perpendicular projection of the field point 
P onto the plane z = 0 . The case shown is where the projection of P lies outside the boundary 
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V, V 



3 

of 2. The “extra” contribution in the sum 2 due to r ^ 2 is cancelled by the final sweep 

i = 1 

in <j> from corner point 4 to comer point 5. Performing the r - integration yields 


H(l, 1,3) = 


~r 

E 

i= 1 



1 



(G.4) 


The next step is to convert equation (G.4) to a line integral along the boundary of Z. To 
do this, consider figure G.2 (which is based on figure D.3). 


Along typical side L of the boundary, r^ = a^ + and the variable of integration is 
related to C as follows: 


cos (f) = 


a 


sin <l> = 


g sign (a) 
+ a^ 


tan (^ = 8 / a 


d0 = 


adC 
+ a^ 


(G.5) 
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Figure G.2. — Geometry Relating to Equation (G.6). 
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Separating equation (G.6) into its two parts gives 




and 



i 



Note: The integral denoted by I2 may be found in reference 28 (page 49, line 3). 
Using these last two results, equation (G.8) becomes 

4 *^i+l 

i = 1 fi. 

where 


P = tan'^ 



tan 


-1 


|h| g \ 
Ia| s/e2 + gV 


In this form, four arc tangents must be computed for each side of 2. 


(G.8) 


(G.9) 

i+' 

(G.IO) 


(G.ll) 

(G.12) 
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To make for more efficient computations, the two arc tangents in equation (G.12) are 
combined into a single arc tangent. 

, |a|E (v/g2 + g2 - |h| ) 

^ = tan-1 ^ 3 y- 

a2 >yc2 + g2 + j2|j,| 


= tan 


-1 


lalE 


g^ + 


Ih| 


(G.13) 


+ g" 


With the aid of the following sketch, 



and the relation 


^ -1 . , -I 

— tan ^ m = tan m 

Im| 


we obtain 




T 

z 


C; 


i+1 


i = 1 ’‘i 


C; 


(G.14) 


(G.15) 
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where 


sin P = 


a £ 

g {j2^ + g2 + lh|) 


COS P = 


+ |h| s/g^ + g^ 

g(s/e2 + g2 + |h|) 


tan /3 = 


a £ 


g^ + Ihl v/c“ + g^ 


(G.16) 


An additional efficiency is gained by combining the difference )3j ^ I - jSj into a single arc 
tangent. For simplicity, consider only a single side of E, with endpoints i = 1 and i + 1 = 2. 

Then 


,r^h\ 

- j3i = tan ^ { ) - tan M — 

\"2/ V^l; 


(G.17) 


where 


C| = g"- + |h| S| 


Sj = + g2 


C-) = g^ + |h| S') 


S2 = 


(G.18) 


Equation (G. 17) reduces to 


i32-j3i = tan' 


ciC2 + a^CiC2 


(G.19) 
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and equation (G.15) becomes 


4 

j= 1 

uL^ 

j = 1 

where the sum is over the four sides of S and the subscripts 1 and 2 are now taken as the 
first and second endpoints of each side j. 

Computing equations (G.19) and (G.20) with the single argument FORTRAN ATAN 
external function returns values of 1^2 ~ in the range (“ 2’ f) ’ double argument 

ATAN2 external function, so as to obtain values of ^2 ~ *n the range (- tt, tt), we must 
also compute siniPi - j3j) and cos (jS-, - )3j ) . These quantities are obtained from equations 
(G. 16), using the difference formulas for sin and cos. The result is 


i^C2Ci -81C2) 
CjCo +a^C|C2 


(G.20) 


sin 



^(^2^1 -^1^2) 

g^djd2 


cos 



c JC 7 + a^CjC2 
g^djd2 


(G.21) 


where 

d| = sj + |h| 

d2 = S2 + |h| (G.22) 


Equation (G.20) is then rewritten as 


4 

H(l,l,3)=j^l ^ tan’’ [sin ^/32 -^i], cos^^ 2 ‘<^i)] 
1 


(G.23) 
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If g"'d|d 2 ^ 0, equation (G.23) can be written as 


H(l, 1,3) = 


|h| 


tan 


-1 




(G.24) 


which is the form appearing in equation (D.41). 

The quantity g h“ is illustrated in figure G.3 for (x, y, 0) ^ 2. Recall from equations 

(D. 1 5) that h = z - Zq where Zq = ax^^ + b and (Xq, y^) is the point on 2 closest to 
( X, y, 0 ). 


P(x,v,z) 



{.X 

2 


Figure G.3. — The Quantity g =y/a^ +h^ for (x,y,o) 4 2 


Figure G.4 shows g =vp+ h^ for the case (x, y, 0) e 2; here, g = 0 if P is a p oint on the 

edge of S, S being defined by f = a^^ + brj^, (^, t?) e 2. [(The quantity s = v'g^ + ^ is also 
shown in figure G.4 since it will be referred to later (in eq. G.33)]. 
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For g = 0, both a and h are zero. Although H( 1 , 1 , 3) is singular for h = 0 , hH( 1 , 1 , 3) is not, 
and it is hH(l, 1,3) that appears in the influence coefficient equations. For the flat panel, 
this can be seen directly from equations (D.l 19), (D.124) and (D.140); for the curved panel 
the corresponding equations are (D.28), (D.30) and (D.33), respectively. An alternate form 
for equations (G.21) can be obtained which handles the g = 0 case provided ^ 
is shown next. 

With considerable algebraic manipulation, the g^ term can be removed from the denominator 
of equations (G.21 ). The result is 
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d|d-> sin 




Ih|(?22 - e,2) 

{h^i *2) 


djd2 cos^/32 ' P]j = + |h| ^s, +$ 2 ) + CjC2 + 


^sjS2 + C|C2^ 


(G.25) 


For (x, y, 0) ^ 2, e.g., as shown in figures G.2 and G.3, the origin of C is on an extension 
of L rather than on L itself. Hence, £j and £2 both of the same sign and £,C2>0. 

For this condition, the terms 82s j + !2jS2 and sjS2 + C1C2 appeanng in the denominator of 
equations (G.25) are never zero (even for g = 0) and djd2 is also positive (d jd2 = Cj£2 for 
g = 0). Thus, equations (G.25) present no computational difficulties for any field point P 
satisfying (x, y, 0) ^ 2. 

To summarize, hH(l, 1, 3) is given by equation (G.23) in terms of the arguments sin ()32 -pj) 
and cos (^2 “ )• For (x, y, 0) e 2 but not on an edge of S (g #0), the arguments are those 

given in equation (G.24). For (x, y, 0) 4 2, equation (G.24) is still valid provided g # 0; for 
this case however, equation (G.23) with arguments given by equations (G.25) is preferred 
since these arguments are valid for any g. Thus, equations (G.23) - (G.25) cover all cases 
except when P is on an edge of S. This case is considered next. 

Consider figure G.5 which shows the field point P approaching Q at an edge of S, along a 
fixed direction given by d. The pertinent geometric quantities g, a, and h are shown for 
two positions of the field point. As P moves to P' to Q along 0 = constant, G moves to 
G' to Q and , 

h - h - 0 

a -*• a' -*• 0 

g - g' - 0 (G.26) 


Now consider equation (G. 1 6) viz., 


tan P 


a£ 

g2+ |h| Jfl- + g2 


(G.27) 
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'i 


f, Z 



Figure G.5. — P Approaching Q Along a Constant Direction d 


As g -»• 0 due to P ->■ Q along B = constant, equation G.27 becomes 


tan j3 


L 1 

|h| |C| 


(G.28) 


From figure G.5, we see that the ratio a/h is equal to tan (y - 0 ) when P reaches Q. 
Equation (G.28) however, contains jy| , so we must distinguish between the case where P 

approaches Q from above, i.e., h > 0, and the case where P approaches Q from below, 
i.e., h <0. These two cases are shown in figure G.6 along with the corresponding angles 
0y and 0g which are both defined over the interval (0, n). 
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With 6 now considered to be either or we can then write 

tan (tt / 2 • ^ , 0 < ^ < tt (Q 

Thus, equation (G.28) becomes 

tan P sign (C) tan (4" (G.30) 

Evaluating this at the two endpoints of L, viz, at C| and C 2 » ^^id noting that C 9 is positive 
and 12] is negative for P at some point on the interior of the edge of S, gives 

Pi ~ P\ 7T- 26 , O<0<7T (G.31) 

where 6 is either 6^ or 6^. 
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Substituting this for one of the sides in equation (G.20) gives 


hH(l. 1,3)= 

|nl 


where the last term is for the side at which P is on the edge of S. 

Equation (G.3 1 ) shows that for P being a point on the edge of S, the value 
edge depends on the direction of approach. Thus, as stated in section D.5, the sum of the 
arc tangent terms making up hH(l, 1, 3) is bounded, but not continuous for P on an edge 
of S (or S for a flat panel); that is, the value of hH(l, 1, 3) is indeterminate. This difficulty 
is avoided for the network edge control points by withdrawing these control points slightly 
from the edge as described in appendix C. 

Finally, we determine the character of hH(l, 1, 3) for the case h = 0, but P is not on an 
edge of S. For (x, y, 0) e 2, this places P at a point on the interior of S. Thus, referring 
to figure G.7, we want to see how hH(l, 1,3) behaves as P approaches the point z^. 

First, write equation (G.16) in terms of the angles a and 0 defined in figure G.7, viz., 

« g 


tan a = 

g 

cos a = — 
s 



sin 9 = 

[hi 

a 

, cos 0 = - 

, 0 < ? < 

(G.33) 


g g 


where 9 is taken as either 0^ or Oj^ in the same fashion as and of figure G.6. So, 
equation (G.16) becomes 


E 


3 sides 


tan 


a(8'>ci - C|C2 

-1 I zJ + (7T-20) 

cjc^ + ^^^1^2) 


(G.32) 


tan P 


a/g 


1+- \/l + (W g)^ 
g 


= (tan a) 


= (tan a) 


cos 6 


1 + sin 0 sj 1 + tan^ a 


^ V 

cos a cos 0 \ 
cos a + sin 0 / 


(G.34) 
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For the case shown in figure G.7, viz., (x, y.O ) = (x^, y^, 0) e S, we see that 6 goes to 
zero as h goes to zero and that a takes on the same value as <p. Hence, equation (G.34) 
becomes 


tan /3 


tan <t> 


h - 0 


(x, y,0) e 2 


(G.35) 


This result is valid for each side of S, so with the aid of figure G.8 we see that equation 
(G.20) yields 


hH(l, 1, 3) I = 

h - 0 
(x,y,0)eS 


|h| 


j = l 




h 
n — 
|h| 


(G.36) 



Figure G.8. — Geometric Interpretation of Equation (G.36) 
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P(x,y,z) 



Figure G.7, — Field Point P Approaching the Point zq on the Interior of S 


173 






If the point P(x, y, z) is located such that (x, y,0 ) ^ 2, e.g., as in figure G.3, the 0’s for the 
four sides of £ do not all go to zero when h goes to zero. This is shown in figure G.9. 

(Note that the lengths ^ are in the plane z = zq whereas the actual a.- are defined in the plane 
z = 0. see figure G.3.) The sid^“facing” the point (x, y, zq) is called side one. We see that 
6 1 is greater than ir/2 and the 0 's for the remaining sides are less than irp. Thus, as h goes 
to zero, 0 1 ->■ 7T and 0; 0 for j = 2, 3, 4. Substituting these values for 0 into equation 

(G.34) gives 


tan /3 


- tan 4> for side I 

+ tan0 for sides 2, 3, 4 


h - 0 


(x, y,0 ) ^ 2 


(G.37) 


P(x,v,z) 



Figure G.9. - The Case h-*0 for (x,y,0) i 2 
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Equation (G.20) then becomes 


hH(l, 1,3) 

h - 0 
(x, y.O) ^ 2 


|h| 




' i‘2 


From figure G. 1 0, we see that 

(^2 ■ *^’l)side ~ ^ (clockwise) 


(G.38) 


^ ^ - 0]^j = r (counterclockwise) 


j = 2 


so 


hH(l, 1, 3) I =0 
h 0 
(x, y,0) ^ 2 


(G.39) 


(G.40) 



Figure G. 10. — Geometric Interpretation of Equations (G.38) — (G.40) 
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Taken together, equations (G.36) and (G.40) give 


(x, y, 0) e £ 
(x, y,0) ^ 2 


hH(l, 1, 3) I 

h - 



(G.41) 


provided P is not on an edge of S (recall equation (G.32)). This result holds for either a 
curved or a flat panel. For a flat panel, S becomes 2, and h = 0 corresponds to the plane 
z = 0 (which contains 2). The effect of the hH(l , 1, 3) jump property (across S) on the 
source velocity, and on the doublet velocity and potential, is discussed in section D.5 for 
the flat panel case. 

G.2 PROCEDURE 1 RECURSIONS: EQUATIONS (D.41 ) ^ (D.48) 

The equations given in procedure 1 are obtained from recombinations of equations (D.37), 
(D.38), (D.39), and (D.56). These four equations are repeated here for convenience 

H(M + 2, N, K) + H(M, N + 2, K) + h2n(M, N, K) = H(M, N, K - 2) (D.37) 


4 

(K - 2) H(M, N, K) = (M - 2) H(M - 2, N, K - 2) - ^ F(M - 1 , N, K - 2), (D.38) 

1 

4 

(K- 2) H(M, N, K) = (N - 2) H(M, N - 2, K - 2) - ^ v F(M, N - 1 , K - 2) (D.39) 

1 • 


p^F(M + 1 , N, K) + v^F(M, N + ] , K) = aF(M, N, K) (D.56) 


Our first task is to put the recursions (D.37), (D.38) and (D.39) into normal form, i.e., 
derive three equivalent recursions in each of which only one index varies. To do this we 
replace M by M + 2 in equation (D.38) and N by N + 2 in equation (D.39) to obtain 

4 

(K - 2) H(M + 2, N, K) = M H(M, N, K - 2) - Z F(M + 1, N, K - 2) (G.42) 

1 

and 

4 

(K - 2) H(M, N + 2, K) = N H(M, N, K - 2) - Z Vj} N + 1, K - 2) (G.43) 
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Substituting H(M + 2, N, K) from (G.42) and H(M, N + 2, K) from (G.43) into equation 
(D.37) and rearranging we obtain 

(K - 2) H(M, N, K) = (K - M - N - 2) H(M, N, K - 2) 

4 4 

+ Z' «^fcF(M+ 1,N, K-2)+ Z i; F(M, N + l,K-2) 
1 1 


Substituting D.56, we obtain a recursion which involves variations in the K index only: 


4 

(K - 2) h2 H(M, N, K) = (K - M - N - 2) H(M, N. K - 2) + Z a F(M, N, K - 2) (G.44) 

1 


Substituting H(M, N, K) from equation (D.38) into (G.44) and rearranging we obtain 
(K - M -N - 2) H(M, N, K - 2) = (M - 2) H(M - 2, N, K - 2) 


-h^ X F(M - 1,N, K -2) - Z aF(M, N, K-2). 
1 1 


Replacing K by K +2 we obtain a recursion which involves variations in the M index only: 


(K - M - N) H(M, N, K) = (M - 2) H(M - 2, N, K) 

4 4 

-h^ Z I'f F(M-1,N, K)- Z iF(M, N, K) (G.45) 
1 1 

Interchanging the roles of M and N, and and we obtain a recursion which involves 
variations in the N index only: 

(K - M - N) H(M, N, K) = (N - 2) h^ H(M, N - 2, K) 


4 4 

-h^ X F(M, N- 1, K)- X aF(M, N, K). (G.46) 

1 ' 1 


Now it is an easy matter to derive equations (D.41 ) through (D.48). Setting M = N = 1 
in equation (G.44) yields equation (D.42). Setting K = 3 in addition and substituting 
H(l, 1, 3) from equation (G.24) yields equation (D.41). Equation (D.43) is obtained by 
setting K = 1 and M = 2 in equation (G.45). Here we note that the first term on the right 
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of (G.45) has a zero coefficient, hence H(0, N, 1 ) is not needed. Equation (D.44) is obtained 
by setting M = 1 and K = 1 in equation (G.46). Equation (D.45) is obtained by setting K = 1 
in equation (G.45). Equation (D.46) is obtained by setting M = 1 in equation (D.39). 
Equation (D.47) is the result of setting M = 2 in equation (D.39). Finally, equation (D.48) 
is simply a rearrangement of the terms of equation (D.37). 

G.3 PROCEDURE 4 RECURSIONS: EQUATIONS (D.61) (D.67) 

The equations given in Procedure 4 ate obtained from recombinations of equations (D.55), 
(D.56) and (D.57). These three equations are repeated here for convenience. 

F(M + 2, N, K) + F(M, N + 2, K) + h^ F(M, N, K) = F(M, N, K - 2) (D.55) 

r-j F(M + 1 , N, K) + F(M, N + 1 , K) = a F(M, N, K) (D.56) 


- (M - 1 ) F(M - 1 , N, K - 2) + (N - 1 ) F(M, N - 1 , K - 2) 

+ (K - 2) F(M + 1, N, K) - (K - 2) F(M, N + 1, K) = E(M, N, K - 2) (D.57) 

Our first task is to put the recursions (D.55), (D.56) and (D.57) into useful form, i.e., derive 
equivalent recursions which will allow systematic evaluation. To do this we solve 
equations (D.5^ and (D.57) for the unknowns F(M + 1, N, K) and F(M, N + 1, K) using 
the fact that = 1 and obtain 

F(M+1,N, K)= [ 2) [(K-2)i'jaF(M, N, K) + (M- j j^_2) 

- (N - 1) F(M, N - 1, K - 2) + E(M, N, K - 2)j (G.47) 

and 

F(M, N + 1, K) = [(K - 2) Ur,a F(M, N, K) - (M - 1) F(M - 1, N, K - 2) 

+ (N- i,K-2)-i'^E(M,N, K-2)J (G.48) 

Replacing M by M + 1 in (G.47), and N by N + 1 in (G.48) we obtain 

F(M + 2, N, K) = [(K - 2) a F(M + 1 , N, K) + M p(^^ N, K - 2) 

-(N - \ F(M + 1, N- 1, K-2)+ E(M + 1, N, K -2^ (G.49) 
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F(M, N + 2, K) = " 2) a F(M, N + 1, K) -(M- 1) v^i^^F(M- 1,N+ l.K-2) 

+ p(m, N, K - 2) - E(M, N + 1 , K - 2)] (G.50) 

Adding equations (G.49) and (G.50) and noting thatv^^ ^ ^ ' 

F(M + 2, N, K) + F(M, N + 2, K) 


= nd^) ' jVjF(M+l,N, K) + i^^F(M, N + l,K)j 


+ (N + M - 1) F(M, N, K-2) 


-(N- Di^tj 
- (M - 1 ) 


Vj F(M + 1,N- l,K-2) + i^^ F(M, N. K - 2) 
F(M, N, K - 2) + F(M - !,N+ l,K-2) 


+ Vjj E(M + 1 , N, K - 2) - E(M, N + 1 , K - 2) 


(G.51) 


The term on the left side of (G.5 1 ) is equal to F(M, N, K - 2) - h F(M, N, K) from equation 
(D.55) . The terms in square brackets on the right side of (G.51) can be reduced via 
equation (D.56). For the first square bracket we use (D.56) directly. For the second, we 
replace N by N - 1 and K by K - 2. For the third, we replace M by M - 1 and K by 
K-2. Making the substitutions and rearranging, we obtain the following recursion on K: 


F(M,N,K)= ^ 

r» ^ 


1 


(K-2) 


- M - N - 


1)F(M, N, K-2) 


+ (N- l)ai;^F(M,N- l,K-2) 


+ (M- l)ai'jF(M- 1,N, K-2)-v^E(M+ 1,N, K-2) 


+ E(M, N + l,K-2)j 

Next we solve equation (D.56) for F(M + 1 , N, K) and for F(M, N + 1 , K): 
F(M + 1.N,K) = ^ F(M,N, K)-^ F(M,N+1,K) 


and 


F(M, N + 1. K) = — F(M, N, K) --p F(M + 1, N, K) 
*'17 V 


(G.52) 


(G.53) 


(G.54) 
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(G.55) 


Replacing M by M + 1 in equation (G.53) we have 

F(M + 2, N, K)=-j^ F(M + 1, N, K)-|^ F(M + 1, N+ 1, K). 

Next we substitute the expression for F(M + 1, N, K) given by (G.53) into the first term 
on the right side of (G.55) and the expression for F(M + 1 , N + 1 , K) given by (G.53) with 
N replaced by N + 1 into the second term on the right side of (G.55) to obtain: 

F(M + 2. N, K) = ^ fa“F(M, N, K) 

- 25 F(M, N + 1 , K) + F(M, N + 2, K)j (G.56) 

Finally, we solve equation (D.55) for F(M + 2, N, K) and substitute the result into the left 
side of (G.56). Upon rearrangement we obtain 

F(M, N + 2, K) = 2 a F(M, N + 1 , K) 

- (a- + F(M, N, K)+ F(M, N, K - 2) (G.57) 

Interchanging the roles of M and N, and and p^, we have the symmetric result; 

F(M + 2, N, K) = 2a p^ F(M + 1 , N, K) 

- (§2 + p^\^) F(M, N, K) + Vr)~ N, K - 2) (G.58) 

We need one more result. Replacing N by N - 1 in equation (G.47) we have 
(K - 2) F(M, N, K) = (K - 2) p^ a F(M, N - 1, K) 

-(M- F(M- 1, N- 1, K-2) 

+ (N - 2) p(]^ N _ 2^ K - 2) 

-i^jE(M, N-l,K-2) (G.59) 

The second and third terms on the right side of (G.59) can be transformed into 
expressions involving K rather than K - 2 using equation (D.55). For the second term, 
we replace M by M - 1 and N by N - 1 in (D.55) and for the third term we replace N 
by N - 2 in (D.55). Upon substitution into (G.59) we obtain: 
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(K-2) F(M, N, K) = (K-2)i;^aF(M, N- l,K)-(M- 1 )j^^j^^F(M+ 1,N- 1,K) 
-(M-l)j'^i'^F(M-l,N + l,K) 

- (M - F(M - 1, N - 1, K) 

+ (N-2)»^j2 F(M + 2, N-2, K) 

+ (N - 2) f(M, N, K) 

+ (N - 2) F(M, N - 2, K) 

N- l,K-2) (G.60) 

Next we replace N by N - 1 in equation (G.53) and substitute the result into the second 
term on the right of (G.60). In addition, we replace M by M - 1 in equation (G.54) and 
substitute the result into the third term on the right of (G.60). Finally, we replace N by 
N - 2 in equation (G.56) and substitute the result into the fifth term on the right of 
(G.60). Making these substitutions and rearranging, we obtain: 

F(M, N, K) = (k_m-N+ 1) [(K - M - 2N + 3) ai^^ F(M, N - 1, K) 

+ (N - 2) (a^ + f(M, N - 2, K) 

-(M- l)ai^^ F(M- 1,N, K) 

- (M - 1 ) F(M - 1 , N - 1 , K) 

-r-^E(M, N- l,K-2)] (G.61) 

Interchanging the roles of M and N, and and we obtain 
F(M, N, K) = (K-M-N + T ) [(K - 2M - N + 3) ai^^ F(M - 1, N, K) 

+ (M - 2) (a- + h^) F(M - 2, N, K) 

-(N- Dai-^ F(M,N- 1,K) 

-(N- F(M- 1,N- 1,K) 

+ i^^E(M- 1,N, K-2)j (G.62) 
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Notice the + sign in the last term on the right side of (G.62) as compared with the - sign 
in the analogous term of (G.6 1 ). This is because the left side of equation (D.57) undergoes a 
change of sign when M and N, and and v are interchanged, hence our formulas will be 
correct if we replace E by - E whenever such an interchange is made 

Now we are prepared to derive equations (D.61) through (D.67). Equation (D.61)is obtained 
from equation (G.52) by setting M = N = 1 . Equation (D.62) is obtained from equation (G.61 ) 
by setting M - K - 1. Equation (D.63) is simply a rearrangement of equation (D.56) with 
K = 1. Equation (D.64) is obtained from equation (G.62) by setting N = K = 1. Equation 
(D.65) is simply a rearrangement of equation (D.56) with K = 1. Equation (D.66) is obtained 
from equation (G.61 ) by setting M = 1 and N = 2. Finally, equation (D.67) is obtained from 
equation (G.57) by setting M = 1 and replacing N by N - 2. 

G.4 PROCEDURE FOR EVALUATING THE E FUNCTIONS 
The E functions are defined by equation (D.58) which we repeat here for convenience: 


E(M, N, K) = 1 1 x)*^" |^(T)-y)^~ ^ 


’ P=V(t-x)^ + (Tj-y)2 + h2 (D.58) 


The E functions required for Procedure 4 can be evaluated recursively using equation 
(D.59) which we also repeat here for convenience: 


P(l)=(x2 + xi) P(I- i)_xjX2 P(I-2) 


(D.59) 


where 



P(I) = A2X2^ - 1 

- A V I- 1 
AiXi 


Examination 

of Procedure 4 indicates that we need to calculate the following E functions: 

a. 

E(2, 1, K-2) 

K = 3, MXFK, 2 


b. 

E(l,2, K-2) 

K = 3, MXFK, 2 


c. 

E(1,N- 1,- 1) 

N = 2, MXQ 


d. 

E(M- 1, 1,-1) 

M = 2, MXQ 


e. 

E(l, l,K-2) 

K = 3, MXK-2, 2 

(G.63) 

These functions can be evaluated using equation (D.59). For this purpose we set E 
Xj , X 2 , Aj , and A^are defined respectively as: 

= P where 
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a. X 2 =P 2 ~"» A 2 = (^2 " /P2’ Aj = (^j - x)/p j ; 1=1, 

b. x-) = P 2 ~“’ A 2 = (i?2 “ y)/P2> A] = Ct7j - y)/p j ; 1=1, 

c. X2 = (t? 2 " y^’ X] = (T?]-y), A2 = P2> Aj=pj; 1=1,MXQ-1 


MXFK - 3 
2 

MXFK - 3 

2 


d. X2 = (£2 " ’ X]=(^j-x), A 2 - P2> Aj-pj; 1-1,MXQ-1 

e. xt = P 2 ~“, A2=1/P2> A| = l/p|; 1=1, ^ (G.64) 

G.5 PROCEDURE 6 RECURSIONS: EQUATIONS (D.99) ^ (D.I02) 

The equations in Procedure 6 are obtained from recombinations of equations (D.96), 

(D.97) and (D.98), which we repeat here for convenience: 

G(M + 1 , N) + G(M, N + 1 ) = a G(M, N) (D.96) 

(D.97) 


(D.98) 


Equation (D. 1 00) is obtained from equation (D.96) by replacing M by M - 1 and rearranging. 
Equation (D. 1 02) is obtained from equation (D.96) by replacing N by N - 1 and rearranging. 
Equation (D.99) is obtained from equation (D.97) by setting M = 1 and replacing N by 
N+ 1. Equation (D. 101) is obtained from equation (D.97) by setting N = 1 and replacing M 
by M + 1 . 

We note from equations (D.99) and (D.lOl) that the D functions must be evaluated for 
the two cases: 

a. D(1,N+1) N=1,MXQ 

b. D(M+1,1) M=1,MXQ (G.65) 


(M- 1)>».^G(M- 1, N) + (N- l)r'^ G(M, N - 1) = D(M, N) 


D(M, N) = * 
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This can be done recursively using equation (D.59). We have 


a. X 2 = t? 2 , A 2 = j? 2 , 1=1, MXQ 


b. X2-^2’ ^l=fl> A2 = ^2> A,=J,; 1=1, MXQ (G.66) 

G.6 CONTINUITY PROPERTIES OF THE H AND F INTEGRALS 

In this section we prove that h^ H(M, N, K) is bounded as a function of P(x, y, z) if 
J = K - M - N and is continuous everywhere and equal to zero when h = 0ifJ> K-M-N. 
In addition, we show that g F(M, N, K) is continuous everywhere and equal to zero when 
g = OifJ>K-M-N+ 1. 

We first note that 


P>lf-x| 

p>|T?-y| 

P > lhl (G.67) 

where 

P = \/(f-x)2 + (T?-y)2 + h2 

It follows that 

p^>|h|^'* 1 ^ P^ where J = K-M-N. (G.68) 

Hence we have 


h-^H(M, N, K) < 




If-x 


1^' ^ |j?-yl^‘ ‘ 


.K 


d ?dt? < !h| 


S P 

= |hH(l, 1,3)1 

From equation (G.24) we see that |h H( 1 , 1 , 3) |< 47 t. 




dfdij 


2 P- 


(G.69) 


Hence 



This shows that h^ H(M, N, K) is bounded when J - M + N - K. Now from equation (D.25) 
it is clear that for arbitrary J, h-^ H(M, N, K) is continuous everywhere except perhaps at 
h = 0. However, we have 

h^ H(M, N, K) = h^ “ [h^ - M - N K)j (G.7 1 ) 

The term in square brackets on the right is bounded by 47 t from the result above, hence, if 
J > K - M - N the expression on the right tends to zero uniformly as h -► 0. Thus 
h^ H(M, N, K) is also continuous at h = 0 and equal to zero there. 

Before considering the F integrals we prove one additional result concerning the H integrals, 
namely that when either M or N is even, then H(M, N, K) is continuous everywhere, 
except when h = 0 and (x, y, 0) belongs to the perimeter of 2. Without loss of generality, 
we assume that M is even. When M = 2, the result follows from equation (D.38) by noting 
from equation (D.40) that all F integrals are continuous everywhere, except perhaps when 
h = 0 and (x, y, 0) belongs to L. For M > 2, the result follows inductively from (D.38). 

Next we show that g:^ F(M, N, K) is continuou s everyw here and equal to zero when g = 0 
ifJ>K-M-N + 1. From the fact that p = + we have in addition to the inequalities 

(G.67), the inequality 

p > g (G.7 la) 

Hence 

pK>|g|J- 1 It-xl'^- 1 |T?-yl^‘ * P ifJ = K-M-N + 2 (G.72) 


It follows that 
IgJ F(M, N, K) 




,J|f- xl^~ H-yl 


N- 1 


X 



= gF(l, 1, 1) 


From equation (D.60) it is easily shown that 

1F(1, 1, 1)1 <Ci +C2lln gl 


(G.73) 


(G.74) 


for some constants C] and C 2 independent of g, hence 

|g^ F(M,N, K)t<g(Ci + C2lln gi)if J = K-M-N + 2 


(G.75) 


This shows that Ig-* F(M, N, K)l tends to zero uniformly as g 0. From equation (D.4U), ^ 

it is c lear that for arbitrary J, g^ F(M, N, K) is continuous everywhere except perhaps at g - 0 where 

p =-\Jg^ + may be zero. However, we have 

gJF(M,N,K) = gJ-(K-^^-N + 2) [gK-M-N + 2p(M,N,K)] 


(G.76) 
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The term in square brackets tends uniformly to zero as g -► 0 from the result above so that 
F(M, N, K) does likewise when J>K-M-N + 2. Hence, g*^ F(M, N, K) is continuous 
everywhere and equal to zero when g = 0 ifJ>K-M-N+l. 


G.7 PROPERTIES OF H* 


In this section we discuss the properties of H* as defined by equations (D.5 1 ) and (D.52), 
which we repeat here for convenience 

H* (M, N, K) = H(M, N, K) - e(M, N, K) (D.5 1 ) 

e(M, N, K) = 27 t KM, N, K) |h|^ + N - K 


if M or N is even 


KM, N, K) = 


[1-1-3-5 • • ■ |M-2ll[M-3-5- ■ •|N-2|] 
[(K-2)(K-4)(K-6) ••• (K-M-N)] 


(D.52) 


First we show that H* satisfies the same recursions as H, in particular that H* satisfles the 
recursions of Procedure 1 for h ^ 0 with the exception of the initial condition (D.41). 
For this purpose, it is sufficient to show that e(M, N, K) satisfies the homogeneous recur- 
sions (D.42) (D.48), i.e., the recursions with all F functions set identically zero. For 

this purpose we note from equation (D.52) that e(M, N, K) satisfies: 


e(M,N,K-2) = 


(K-2)e(M, N, K) 
(K-M-N-2) 


(G.77) 


c(M -I- 2, N, K) = 


h^ Me(M, N, K) 
(K-M-N-2) 


(G.78) 


e(M,N-t-2, K) = 


h^ Ne(M, N, K) 
(K-M-N-2) 


(G.79) 


We note that the initial condition 


6(1, 1, 1) = 27t| hi (G.80) 

Then equation (D.42) (with F = Oand H replaced by e) follows from equation (G.77) with 
M = N = I . Equation (D.43) follows from the (D.52) and the fact that M is even. Equation 
(D.44) follows from equation (G.79) with M = K = 1 and N replaced by N - 2. Equation 
(D.45) follows from equation (G.78) with K = 1 and M replaced by M - 2. Equation (D.46) 
follows from equation (G.79) with M = 1 and N replaced by N - 2 combined with equation 
(G.77) with M = 1 and N replaced by N - 2. Equation (D.47) follows from 
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equation (D.52) and the fact that M is even. Finally, equation (D.48) may be obtained by 
subtracting (G.77) from the sum of (G.78) and (G.79), and then replacing M by M - 2. 

Next we prove that H*(M, N, K) is continuous when (x, y, 0) belongs to the interior of 2. 
From section (G.6) we know that H(M, N, K) is continuous everywhere when K < M + N 
and bounded when K = M + N. However, when K = M + N, either M or N must be even since 
K is odd, hence, we know in addition that for K = M + N, H(M, N, K) is continuous every- 
where except perhaps when (x, y, 0) belongs to the perimeter of 2 and h = 0. It follows 
from equation (D.5 1 ) that H*(M, N, K) is continuous when (x, y, 0) is in the interior of 2 
and K < M + N. Hence, we consider only the case K > M + N. For this case we shall show 
that 


K) = - JJ^ ^i- 


x)^ ~ ^ (t? - y)^ “ * 


dfdi?. 


(G.81) 


where Z * is the exterior of S, i.e., the whole tj plane minus the quadrilateral Z. From 
equations (D.51) and (D.35) it is sufficient to show that 


e(M,N,K)= JJ" 

J - T 7 Plane 






dfdT7 


(G.82) 


To do this we let ({ - x) = |h|r cos0 and (t? - y) = |h|r sin0. Then 
JJ' (f-x)^-^(T?-y)N- 1 
5 - T} Plane 


d^dr? 


= Ih| 


M 


^27T 


+ N-K ^ r (cos«)^!-l(slnO)N-l 

•'O •{) (1 +r2) 


drdd (G.83) 


Hence, it is sufficient to show that 

..27T 


KM, N, K) = 


X 

27T 


r .M + N-i ^ 

j I 

*'o •'o (i + n 


^ drd0 




(G.84) 
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By symmetry, the second integral in square brackets is zero unless M and N are both odd. 
Hence, we need only consider the case when M and N are both odd. 


Itt 



lir 

(cosd) 


M- 1 


(sin0)'^ ' • 




7t/2 

(cosfl)^ “ * 


(sin0)’^ " • 


dO 


(G.85) 


The first integral in square brackets can be put in similar form by setting r = tan 0 to obtain 

.■nil 

J- XT _ 1 1/ X/f XT 1 

(G.86) 


/ J.M + N - 1 r* 

^ 2^ K/2 ^ J + N - 1 (cos0)^ ■ ^ " N - I 


0 

From integral tables (ref. 29) we have that 
.7 t/2 


f (cos x)^ ' 1 (sin x)^' " ^ dx ~ ^ ^ ^ ^ ^ ^ ^ 


(G.87) 


where F is the F function and 


ii)- 


1-3-5 • • • IM-21 ^ . . 

: ^ — )]2 If Mis odd 


(G.88) 


Applying (G.87) to (G.86) and (G.85), and substituting into (G.84), we need only to 
show that 


KM,N, K) = ^ 


F 1 

1^1 

1 r 1 

© 

r 1 

fK-M-N\ 
2 / 

F 

(!) 

1 


(G.89) 


Upon applying (G.88) to (G.89) and comparing with equation (D.52) we obtain the 
desired result. 


G.8 VALIDITY OF REVERSE RECURSIONS (PROCEDURES 2 AND 5) 

In this section we justify use of the initial condition (D.49) in the reverse recursion (D.50). 
The argument is precisely the same for the use of equation (D.68) as an initial condition 
for equation (D.69) and we therefore consider only equations (D.49) and (D.50) which 
we repeat here for convenience 
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H(l, 1 , NHK + MXK) = 0. 


(NHK= 16) 


(D.49) 


H(l, l,K-2) = 


1_ 

(K-4) 



(K-2)H(1, 1, K) 


2 aF(l, l,K-2) 
1 


(D.50) 


for K = NHK + MXK, 3, - 2. 


Without loss of generality, we assume (x, y, 0) ^ 2. (If (x, y, 0) e 2, we replace H by 
H* in which case we are dealing with (D.49) and (D.50) as applied to Procedure 3.) In 
either case, we have p >d where dj^ is defined as the minimum distance from (x, y, 0) 
to the perimeter of 2. Then 


H(l, 1,K + J) = 


rr di;dv 

jj 

2 ^ 




d^dn _ J_ 
7^ dH^ 


H(l, 1,K) 


It follows from the assumption of Procedure 2 (|h| < dpj) that 
Ih-* H(l, 1,K + J)K||3^ 

Setting K = MXK and J = NHK we have 


H(l, 1,K)|=57 IH(1, 1,K)1 


(G.90) 


(G.91) 


IhNHK H( 1 , 1 , NXK + MXK) | < 6 |H( 1 , 1 , MXK)| 


(G.92) 


By substituting equation (D.49) into equation (D.50) for K = NHK + MXK and then 
successively solving for H( 1 , 1 , K) with lower values of K, we find 


H( 1 , 1 , MXK) = h^HK^d, 1, NHK + MXK) + F terms (G.93) 

Since 6 = .01 we see from (G.92) that the first term on the right of (G.93) is negligible 
compared with H(l, 1, MXK), hence, it can be ignored. This is easily accomplished by 
setting H( 1 , 1 , NHK + MXK) = 0. 
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